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ABSTRACT 

This  Semiannual  Technical  Summary  describes  the  Lincoln 
Laboratory  Vela  Uniform  program  for  the  period  1  April 
through  30  September  1981.  During  this  period,  the  first 
working  prototype  of  a  Seismic  Data  Center  has  been  com¬ 
pleted.  in  this  report,  Sec.  I  describes  this  prototype  sys¬ 
tem,  Sec.  II  describes  a  series  of  activities  in  seismic 
processing  related  to  the  Center,  and  Sec.  Ill  describes  a 
series  of  investigations  in  General  Seismological  Research. 
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SUMMARY 


This  is  the  thirty-fifth  Semiannual  Technical  Summary  (SATS)  report  describing  the 
activities  of  Lincoln  Laboratory  funded  under  Project  Vela  Uniform.  This  report  covers  the 
period  1  April  through  31  September  1981.  Project  Vela  Uniform  is  a  program  of  research 
into  the  discrimination  between  earthquakes  and  nuclear  explosions  by  seismic  means.  An 
important  recent  emphasis  of  the  project  is  in  the  development  of  data -handling  and  analysis 
techniques  that  may  be  appropriate  for  the  monitoring  of  a  potential  C  omprehensive  Test  Ban 
Treaty.  The  FY  1981  program  consists  of  the  development  of  a  Prototype  Seismic  Data  Center 
(SDC),  which  will  fulfill  U.  S.  obligations  that  may  be  incurred  under  such  a  Treaty,  and  which 
will  also  provide  a  data  resource  to  the  seismological  research  community.  The  program  also 
contains  an  element  of  seismic  research,  with  emphasis  on  those  areas  directly  related  to  the 
operations  of  the  SDC. 

During  this  period,  the  first  working  Prototype  SDC  has  been  completed.  This  system  is 
described  in  Sec.  I  of  this  report.  It  is  capable  of  accepting  input  of  both  parameter  and  con¬ 
tinuous  waveform  data,  carrying  out  automatic  signal  detection  and  parameter  extraction  on 
the  waveform  data,  merging  all  parameters,  and  carrying  out  automatic  association  to  produce 
an  event  list.  Review  by  a  human  analyst  is  provided  for,  both  after  automatic  signal  processing 
and  after  automatic  association.  Analyst  interaction  utilizes  a  single-user  interactive  computer 
graphics  system  called  the  Seismic  Analysis  Station  (SAS).  This  Prototype  will  form  the  basis 
for  future  SDC  development. 

In  Sec.  II,  several  studies  in  the  area  of  Seismic  Processing  are  described.  One  study 
demonstrates  that  augmented  transition  network  (ATN)  grammars  can  be  of  significant  use  in 
the  recognition  of  glitch  noise.  Another  compares  ATN  grammars  with  dynamic  time  warping 
as  means  for  the  recognition  and  characterization  of  seismic  signals.  A  third  study  describes 
the  affinity  technique  for  the  analysis  of  the  morphological  structure  and  information  content 
of  seismic  signals.  We  are  beginning  to  apply  the  concepts  of  expert  systems  to  the  problem 
of  automatic  association  (AA),  and  expect  that  this  will  lead  to  much  more  intelligent  AA  pro¬ 
grams.  Another  investigation  examines  the  use  of  hash  functions,  and  the  application  of  en¬ 
tropy  concepts  to  measure  their  performance.  These  functions  may  have  an  important  appli¬ 
cation  to  the  retrieval  of  information  from  database  storage.  We  have  begun  to  study  potential 
improvements  in  signal  detection.  It  is  shown  that  the  number  of  significant  detections  ob¬ 
tained  in  a  given  time  period  can  be  increased  by  carrying  out  signal  detection  on  the  hori¬ 
zontal  components  of  3-component  instruments.  And  we  have  completed  an  extensive  compari¬ 
son  of  the  performance  of  automatic  techniques  and  the  human  analyst  in  picking  the  first 
arrivals  of  seismic  signals.  While  the  automatic  methods  produce  more  false  alarms,  they 
are  comparable  to  the  analyst  in  detecting  real  signals.  Analysts  were  shown  to  be  less 
reliable  than  commonly  assumed. 

Section  III  contains  reports  on  a  number  of  areas  of  general  seismological  research.  A 
detailed  review  of  frequency-dependent  anelasticity  has  been  carried  out.  Additional  evidence 
for  absorption-band  models  is  given.  Another  study  approaches  the  problem  "f  including 
transverse  isotropy  in  ray  theoretic  calculations.  Several  new  extensions  of  a  method,  devel¬ 
oped  earlier,  for  the  simultaneous  estimation  of  seismic  moment  tensor  and  event  hypocenter 
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are  described.  Application  of  these  methods  to  several  recent  earthquakes  is  also  included. 

In  one  case,  evidence  is  given  that  a  complex  series  of  motions  occurred  at  the  seismic  source. 
Finally,  an  application  of  a  new  method  for  the  analysis  of  mantle  waves  shows  how  these  can 
be  used  to  extract  direct  information  on  Earth  structure. 

M.  A.  Chinnery 
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SEISMIC  DISCRIMINATION 


I.  PROTOTYPE  SEISMIC  DATA  CENTER 

A.  SYSTEM  FUNCTIONS 

The  Prototype  Seismic  Data  Center  (SDC)  described  in  this  and  later  sections  of  this  report 

has  been  assembled  at  Lincoln  Laboratory.  The  Prototype  SDC  has  been  based  on  a  Lincoln 

Laboratory  overall  system  design,1  and  various  stages  in  the  development  of  the  system  have 
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been  documented  in  previous  SATS.  The  intent  of  this  Prototype  SDC  is  that  it  should  be  a 
■working  system,  capable  of  carrying  out  all  the  basic  functions  that  will  be  required  in  a  full- 
scale  SDC,  and  able  to  be  used  to  test  out  a  variety  of  design  concepts  and  operational  proce¬ 
dures.  The  system  described  in  this  report  has  been  successfully  demonstrated  to  DARPA, 
and,  with  minor  improvements,  will  form  the  basis  for  further  SDC  development. 

The  functions  that  the  Prototype  SDC  will  perform  are: 

(1)  Accept  input  data  in  the  form  of  both  lists  of  parameters  describing 
seismic  arrivals  and  also  continuous  waveform  data.  Neither  of  these 
data  types  is  currently  available  on-line,  so  the  Prototype  SDC  is  based 
on  input  from  magnetic  tape. 

(2)  Entry  of  both  parameter  and  waveform  data  into  a  suitable  database 
management  system.  The  INGRES  system  designed  by  Lawrence 
Berkeley  Laboratory  has  been  chosen  for  this  task. 

(3)  Automatic  processing  of  the  continuous  waveform  data,  to  both  detect 
seismic  signals  and  also  extract  parameters  describing  those  signals. 

The  outputs  of  this  stage  are  a  list  of  arrival  parameters  and  a  set  of 
waveform  segments  each  containing  at  least  one  detected  signal. 

(4)  Analyst  review  of  the  output  of  the  automatic  processing  algorithms,  in¬ 
cluding  the  elimination  of  false  alarms,  the  identification  of  seismic 
phases,  and  the  refinement  of  the  various  seismic  parameters  mea¬ 
sured,  such  as  origin  time,  amplitude,  etc. 

(5)  Merging  of  the  refined  arrival  data  extracted  from  the  waveforms  with 
all  other  arrival  data  in  the  INGRES  database  management  system. 

(6)  Automatic  association  of  the  arrival  data  into  a  Preliminary- Event 
List  (PEL). 

(7)  Analyst  review  of  the  PEL,  including  changing  arrival  associations, 
adding  depth-phase  information  and  refinement  of  calculated  hypocenters. 

Also,  review  of  those  arrivals  not  associated  with  events  in  the  PEL, 
with  a  view  to  adding  new  events. 

(8)  Preparation  and  publication  of  a  completed  Final-Event  List  (FEL). 

In  the  sections  that  follow,  we  describe  how  each  of  these  functions  is  accomplished. 

M.  A.  Chinnery 
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B.  HARDWARE  CONFIGURATION 

The  hardware  configuration  of  the  Prototype  differs  in  scale  from  the  projected  SDC  and  in 
a  few  details  in  the  subsystems  represented.  This  section  details  the  configuration  in  use  for 
the  Prototype  at  this  time. 

The  logical  configuration,  as  shown  in  Fig.  I-i,  includes  the  Seismic  Analysis  Station  (SAS), 
Database  Subsystem  (DBS),  and  Automatic  Processing  Subsystem.  The  Communications  Inter¬ 
face  Subsystem  is  represented  only  by  routines  used  to  enter  parameter  and  waveform  data  into 
the  DBS.  A  change  from  the  SDC  design  is  the  use  of  a  direct  link  between  the  VAX  11/780  which 
supports  all  the  subsystems  except  the  SAS  and  the  PDP- 11/44  which  implements  the  SAS.  This 
direct  link  uses  the  DEC  DMC-11  high-speed  communication  link  to  send  data  between  the  ma¬ 
chines.  The  DMC  is  used  with  a  command  which  transmits  files  between  the  two  machines.  This 
command,  as  well  as  the  details  of  the  capabilities  of  the  other  subsystems,  is  covered  in  sub¬ 
sequent  sections. 

The  hardware  configuration  is  covered  in  detail  by  Fig.  1-2,  a  schematic  presentation  of  the 
hardware  configuration  used  by  the  Prototype  system.  The  hardware  configuration  is  listed  by 
computer  in  Table  1-1  which  gives  the  part  number  and  supplier  for  all  major  components  of 
both  computers.  A.  G.  Gann 

C.  DATABASE  MANAGEMENT 

This  report  describes  the  database  for  the  SDC  Prototype  system.  The  database  contains 
both  wamform  and  parametric  data  and  is  based  on  INGRES.  INGRES  is  a  relational  database 
management  system.  A  relation  can  be  thought  of  as  a  table  of  information.  In  INGRES  termi¬ 
nology,  a  "tuple"  is  a  row  of  a  table  (relation)  and  a  "domain"  is  a  column.  Commands  to  INGRES 
are  given  in  a  language  called  "Quel."  Quel  commands  can  be  used  to  extract  from,  modify,  and 
do  computations  on  a  relation  or  combinations  of  relations.  INGRES  can  be  accessed  directly 
through  the  INGRES  monitor  or  indirectly  from  a  C  program  using  Equel,  a  C  preprocessor. 
Most  of  the  primitive  database  commands  are  implemented  as  Quel  macros.  Only  the  waveform 
primitives  and  the  primitives  that  perform  seismic  calculations  need  to  be  implemented  in  Equel. 

1.  Database  Structure 

Parametric  data  are  measurements  made  from  seismograms,  such  as  arrival  times,  ampli¬ 
tude,  and  periods.  Such  data  are  determined  by  the  SDC  itself  and  also  from  data  sent  to  it  from 
other  sources.  These  data  will  be  organized  into  events  using  automatic  and  interactive  methods. 

Parameter  data  related  to  an  event  are  organized  into  a  4-level  structure  as  shown  in 
Fig.  1-3.  An  event  consists  of  a  number  of  origins  —  each  origin  being  an  estimate  of  the  event's 
location  based  upon  evidence  from  the  arrivals  with  which  the  origins  are  associated.  One  ori¬ 
gin,  referred  to  as  the  preferred  origin,  will  be  used  to  represent  the  best  estimate  of  the  origin 
of  the  event.  Allowing  an  event  to  have  multiple  origins  is  useful  for  several  reasons: 

(a)  Historically,  data  centers  have  reprinted  epicenter  estimates  provided  by 
other  agencies. 

(b)  Multiple  origins  allow  an  analyst  to  have  several  competing  hypotheses 
for  one  event  during  the  association  process.  Normally,  after  association 
is  complete,  one  origin  is  chosen  to  represent  the  event  while  others  are 
discarded. 
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TABLE  1-1 

DEMONSTRATION  HARDWARE  COMPONENTS 

VAX  1 1/780  (Model  SV-AXHHB-CA),  (4  boys),  with 

1 

LA-120  console 

(DEC) 

4.0-Mbyte  Memory 

(2.5  MB  DEC, 

1.5  MB  Motorola) 

1 

H7U2-A  Battery  backup  for  memory 

(DEC) 

1 

FP780-AA  Floating-point  accelerator 

(DEC) 

2 

DW780-AA  Unibus  adopter 

(DEC) 

1 

BAU-KE  Unibus  extension 

(DEC) 

2 

RK07  28-Mbyte  disk  drives  with  Unibus  controller 

(DEC) 

32-port  terminal  controller,  DZ1 1 

(DEC) 

1 

DMC-lt  1-Mbaud  intercomputer  link  to  SAS 

(DEC) 

1 

600-1 pm  line  printer 

(Southern  Systems) 

2 

STC  tape  drives  (1600/6250  bpi,  75  ips) 
with  AVIV  controller 

(AVIV/STC) 

4 

CDC  675-Mbyte  (nonremovable)  disk  drives 
with  SI  controller 

(SI/CDC) 

2 

CDC  300-Mbyte  disk  pack  drives 
with  SI  controller 

(SI/CDC) 

Dial-in  modem  (300/1200  bd) 

(Rocal  Vadic) 

Datagraphix  (132-column)  video  display  terminal 
for  SAS 

(GD) 

32-V  BSD  4.0  UNIX  operating  system 

(WECo/UCB) 

PDP-11/44  (Seismic  Analysis  Station)  with 

1 

LA-120  console 

(DEC) 

512-kbyte  memory 

(DEC) 

Battery  backup  for  memory 

(DEC) 

1 

FP11-F  Floating-point  processor 

(DEC) 

1 

DDI1-0K  Unibus  expansion 

(DEC) 

2 

RL02  10-Mbyte  disk  cartridge  drives 
with  controller 

(DEC) 

16-port  terminal  controller,  DZ11 

(DEC) 

1 

DMC-11  1-Mbaud  intercomputer  link 

(DEC) 

1 

LA180  line  printer  with  parallel  interface 

(DEC) 

2 

Cipher  tape  drives  (800/1600  bpi,  45  ips) 
with  AVIV  controller 

(AVIV/Cipher) 

1 

CDC  300-Mbyte  disk  pack  drives 
with  SI  controller 

(SI/CDC) 

UNIX  V7  operating  system 

Megatek  7000  refresh  graphics  display  system  with 

M702  controller,  IPCU,  data  tablet,  joystick, 

HCRST  2-dimensional  hardware,  10-mil  monitor, 
extra  memory,  rasterizer  and  interface  for  Versatec 
hardcopy  unit 

Versatec  printer/plotter  (200-dot/in.  resolution) 

(WECo) 

Datagraphix  (132-column)  video  display  terminal 

(GD) 
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(c)  As  new  processing  methods  are  developed,  new  results  can  be  added  to 
the  database  without  erasing  older  published  results. 

Each  origin  contains  information  describing  the  hypocenter  as  well  as  a  list  of  arrivals 
which  are  associated  with  that  origin. 

A  Related  Arrival  Group  (RAG)  is  a  group  of  arrivals  recorded  by  one  station  which  are 
related  to  one  event,  such  as  a  p-  and  an  s-arrival. 

Only  an  index  to  the  waveform  data  is  maintained  in  the  database.  The  actual  waveforms 
exist  as  UNIX  files  in  gram-file  format.  This  allows  some  flexibility  in  managing  disk  space. 
The  connection  between  arrivals  and  waveform  data  is  not  maintained  explicitly  in  the  database. 
Only  when  portions  of  the  database  are  extracted  for  analysis  is  the  correspondence  set  up. 

2.  Database  Primitives 

The  database  primitives  are  only  an  initial  set,  and  are  designed  to  maximize  ease  of 
implementation  at  the  expense  of  ease  of  use.  Once  these  primitives  are  working  satisfacto¬ 
rily,  a  more  convenient  user  interface  can  be  developed. 

To  use  the  primitives,  the  user  establishes  a  working  context  within  the  database  structure: 

EVENTS 

ORIGINS 

RAGS 

ARRIVALS 

One  can  define  the  context  at  any  of  the  four  levels,  with  or  without  respect  to  the  other  levels. 
For  example,  to  work  on  origin  51  of  event  5131,  specify: 

work  event  5131 
work  origin  51 

All  subsequent  commands  will  now  be  with  respect  to  this  context. 

The  analyst  can  define  a  new  RAG  without  first  defining  the  event  or  origin  context,  which 
may  be  unknown: 

work  event  0 

new  rag  (range  =  *1") 

new  arrival  (day  =  80105,  time  =  1013.1,  \ 
phase  =  "pg”,  qual  =  "i") 

new  arrival  (day  =  80105,  time  =  1045. 0,\ 
phase  =  "lgn,  qual  =  "e") 

new  rag  ( ) 

work  rag  0 
work  event  0 

A  primitive  continued  on  another  line  requires  a  backslash,  (\). 
a.  Arrival  Primitives 


new  arrival  (  attribute  =  value,  . . .  ) 

Creates  a  new  arrival  and  makes  it  the  working  arrival 


r 


update  arrival  arrival- id  (  attribute  =  value,  ...  ) 

Modifies  an  existing  arrival 

remove  arrival  arrival-id 

Removes  an  arrival  from  the  database 

add  arrival  arrival-id 

Adds  an  arrival  to  the  working  RAG 

sub  arrival  arrival-id 

Removes  an  arrival  from  the  working  RAG 

pr  arrival  [arrival-id] 

Prints  an  arrival 

b.  RAGs  Primitives 

new  rag  (  attribute  =  value,  . .  .  ) 

Creates  a  new  RAG  and  makes  it  the  working  RAG 

work  rag  rag-id 

Sets  up  a  new  working  RAG 

update  rag  rag-id  (  attribute  =  value,  . . .  ) 

Modifies  an  existing  RAG 

remove  rag  rag-id 

Removes  a  RAG  from  the  database 

pr  rag  [rag-id] 

Prints  information  about  a  RAG 

c.  Origin  Primitives 

new  origin  (  attribute  =  value,  ...  ) 

Creates  a  new  origin  and  makes  it  the  working  origin 

work  origin  origin-id 

Sets  up  a  new  working  origin 

update  origin  (  attribute  =  value,  . . .  ) 

Modifies  an  existing  origin 

copy  origin 

Makes  a  copy  of  the  working  origin 

remove  origin  [origin-id] 

Removes  an  origin  and  its  associated  arrivals 

prefer  origin  [origin-id] 

Identifies  an  origin  as  the  preferred  origin 

pr  origin 

Prints  a  summary  of  the  working  origin 
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d.  Association  Primitives 

assoc  arrival  [arrival- id] 

Associate  an  arrival  to  the  working  origin 

assoc  rag  [rag-id] 

Associate  the  arrivals  from  a  RAG  to  the  working  origin 

unassoc  arrival  arrival-id 

Unassociate  an  arrival  from  the  working  origin 

unassoc  rag  rag-id 

Unassociate  the  arrivals  of  a  RAG  from  the  working  origin 
def  arrival- id 

Makes  the  arrival  identified  by  arrival-id  defining  for  the  working  origin 
of  the  working  event 

nondef  arrival-id 

Makes  an  arrival  nondefining  in  the  working  origin 
nondef  resid  tolerance 

Any  arrival  associated  to  the  working  origin  that  has  an  absolute  residual 
less  than  tolerance  is  made  nondefining 

e.  Event  Primitives 

new  event  (  attribute  =  value. . .  ) 

Create  a  new  event  and  make  it  the  working  event 

work  event  event-id 

Make  an  event  the  working  origin 

remove  event  event-id 

Remove  an  event  and  all  of  its  origin  information 
pr  event 

Prints  a  summary  of  the  working  event  including  the  preferred  origin 
and  arrival  information 

print  olist 

Prints  a  summary  of  the  list  of  origins  for  the  working  event 
evbul  day!  ttmel  day2  time2 

Print  an  event  bulletin  for  the  corresponding  time  interval 
loc  [— o] 

Adjust  the  working  event  location  using  the  working  origin  as  the 
starting  location 

match  (-t  toler]  [-1  interval]  [-p]  [-r] 

Finds  all  arrivals  which  match  the  working  origin  of  the  working  event 
to  within  a  given  tolerance 

f.  Waveform  Primitives 
loadgi 

Adds  waveforms  from  a  gram  database  into  the  waveform  database 
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gdbs  [-IXvsc]  ingresdb  gdbdir  gdbname  start  interval  station  [channel]. . 
Create  a  gram  file  by  station  for  a  given  time  interval 


gdbp  [- IXGPvaoepd]  ingresdb  dir  gdbname  [-s,  +s]  [date]  [id...] 

Create  a  new  gram- file  database  for  arrivals  extracted  by  event,  or 
origin 


K.R.  Anderson 
P.  Kreps 
K.  J.  Schroder 


D.  AUTOMATIC  SIGNAL  PROCESSING 

In  the  Prototype  SDC,  input  continuous  waveform  data  are  subjected  to  a  series  of  automatic 
processing  steps  which  are  summarized  as  follows: 

(1)  Signal  detection,  and  the  extraction  of  a  waveform  segment  containing 
the  detected  signal 

(2)  Estimation  of  signal  onset  time 

(3)  Discrimination  between  local  and  teleseismic  signals,  based  on  the 
frequency  content  of  the  signal 

(4)  Measurement  of  signal-to-noise  ratio  (SNR) 

(5)  Measurement  of  the  direction  of  first  motion  for  impulsive  signals 

(6)  Measurement  of  signal  amplitude  and  dominant  period. 

In  the  first  version  of  the  Prototype  SDC,  input  waveform  data  are  stored  on  disk  on  the 
VAX  (see  Secs.  A  and  B  above)  and  then  retrieved  for  analysis.  The  output  waveforms  segments 
and  signal  parameters  are  then  entered  into  the  DBS.  Later  versions  of  the  Prototype  SDC  will 
utilize  the  DBS  for  both  input  and  output  data.  The  various  automatic  processing  steps  are  de¬ 
scribed  individually  below. 

Signal  detection  is  accomplished  using  a  "Z-statistic"  algorithm  based  on  that  formulated 
5  4 

by  Swindell  and  Snell.  This  detector  has  been  described  earlier.  The  various  parameters 
describing  this  detector  are  listed  in  Table  1-2,  together  with  typical  values  of  these  parameters 
which  are  loaded  at  run  time. 

Following  a  detection,  a  waveform  segment  extending  from  i  min.  before  the  detector  "on* 

time  to  1  min.  after  the  detector  "off"  time  is  passed  to  a  signal  analysis  program.  This  pro- 
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gram  is  based  on  the  algorithm  of  Allen  and  has  been  described  earlier.  This  program  first 
estimates  the  signal  onset  time  by  refiltering  the  input  waveform,  and  then  applying  a  com¬ 
pressed  LTA  (10  s)/STA  (0.1  s)  method.  It  also  measures  signal  power  in  two  frequency  bands, 
and  makes  a  qualitative  estimate  of  event  distance  (local  or  teleseismic).  Current  frequency 
bands  have  been  based  on  SRO  data,  and  may  need  to  be  revised  for  RSTN  waveforms.  For 
teleseismic  events,  amplitude  is  measured  as  the  maximum  observed  in  the  5  s  following  signal 
onset  and  converted  to  millimicrons  of  ground  motion,  and  a  dominant  period  is  measured  for 
this  maximum  peak.  The  quality  of  the  arrival  (impulsive  or  emergent)  is  determined  by  the 
SNR  in  the  vicinity  of  the  onset  time  (if  the  STA/LTA  ratio  is  greater  than  10,  the  signal  is 
designated  impulsive).  When  the  signal  is  impulsive,  a  direction  of  first  motion  is  measured. 

The  signal  analysis  program  then  outputs  a  "new- arrival"  primitive  for  entry  into  the 
INGRES  database.  This  primitive  contains  station,  channel,  date,  onset  time,  quality,  phase 
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TABLE  1-2 

SIGNAL  DETECTOR  PARAMETERS 

Parameter 

Typical  Value 

Bandpass  filter,  low  corner  frequency 

0.9  Hz 

Bandpass  filter,  high  corner  frequency 

3.5  Hz 

Length  of  LTA  window 

45  s 

Length  of  STA  window 

3  s 

LTS  to  STA  lag 

15s 

Interval  between  successive  STA  values 

1.5  s 

Freeze  LTA  between  detections  (yes/no) 

No 

Maximum  LTA  freeze  time 

5  min. 

Initializing  mean  of  filtered  trace 

12 

Initializing  standard  deviation  of  filtered  trace 

0.35 

Detection  threshold  (standard  deviations) 

4 

Successive  triggers  needed  to  define  detection 

2 

Time  after  last  trigger  before  detection  declared  over 

20  s 

identification  (normally  P)  first-motion  direction,  distance  range,  amplitude,  and  period.  An 
author  entry  is  also  included  to  distinguish  parameters  determined  automatically  from  those 
measured  by  an  analyst. 

Future  enhancements  of  automatic  signal  processing  include  the  addition  of  synthetic  anal¬ 
ysis  using  augmented  transition  network  grammars.  With  these  tools  it  will  be  possible  to 
identify  and  remove  glitches,  calibration  pulses,  and  various  kinds  of  nonseismic  noise.  They 
should  also  be  applicable  to  the  problem  of  identifying  S  signals,  particularly  for  local  events. 

M.  A.  Tiberio 
M.  W.  Shields 


E.  INTERCOMPUTER  COMMUNICATION 

Two  major  efforts  have  been  completed  for  data  communication  between  the  VAX  11/780 
database  system  and  PDP-ll/44  of  the  SAS.  The  major  projects  were  installation  of  the  Bell 
Laboratory  UUCP  package  and  development  of  wideband  data  bus  (i- Mbaud  Digital  Equipment 
Corporation  DMR-il)  with  file  transfer  protocol  (Fig.  1-4).  This  hardware  eventually  will  be 
replaced  by  the  local  network  hardware. 

1.  Installation  of  UUCP  Package 

The  purpose  of  UUCP  installation  between  VAX  11/780  and  PDP-11 /44  is  to  provide  low- 
speed  file  copy  and  command  execution.  Gradually,  these  features  will  be  replaced  by  locally 
developed  software.  This  provides  an  interim  communication  mechanism  during  DMR  develop¬ 
ment,  and  provides  some  features  not  yet  implemented  in  the  DMR  system. 
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2.  1  -  Mbaud  DMR-11 

Early  development  used  the  DEC  DMC-11  hardware  because  the  DMR-11  was  not  yet  avail¬ 
able.  A  number  of  hardware  changes  were  made  to  fix  DMC  hardware  problems:  long  DMA 
requests  would  hang  the  Unibus  and  the  DMC  spuriously  issued  RDYO  interrupts  when  running 
full-duplex  under  heavy  loads.  Also,  the  sequencing  and  timing  of  commands  to  initialize  the 
DMC  were  found  to  be  critical.  DEC  claims  that  these  problems  have  been  corrected  on  the 
DMR-11  (Fig.  1-5). 

The  DMR  microcomputer  provides  a  reliable  link,  performing  most  of  the  low-level  pro¬ 
tocol  work  such  as  error  checking  and  retransmission.  The  host  device  drivers  must  provide 
general  consistency  checking  and  flow  control. 

Drivers  have  been  written  for  the  PDP-11  and  VAX  11  Unibus  interfaces  using  DMA  data 
transfer  in  512-byte  packets.  The  latest  tests  show  300-kbaud  bandwidth  process  to  process 
between  two  machines  separated  by  a  DMC-11  network  link.  As  soon  as  DMR  installation  is 
finished,  we  will  try  increasing  the  packet  size  to  achieve  better  performance. 


3.  Easy  File  Transfer  Protocol 

Reliability,  compatibility,  and  ease  of  use  are  the  main  concerns  of  this  protocol.  All 
transfer  requests  are  initiated  by  a  control  program  running  on  the  VAX.  File  transfers 
(Fig.  1-6)  are  then  conducted  by  the  control  program  and  a  server  program  running  at  the  SAS 
via  the  DMR.  The  command  syntax  is  made  to  be  compatible  with  UUCP  and  CP,  except  that 
the  file  name  is  prefixed  with  host  name.  The  best  Unix  file-to-file  transfer  rates  are  about 
50  kbps. 

The  server  program  running  at  the  Analyst  Station  performs  like  a  function  dispatcher. 

It  simplifies  adding  future  functions,  e.g.,  virtual  file,  command  execution,  to  the  structure. 

H.  Hsiung 

F.  WAVEFORM  DISPLAY 

Work  has  progressed  on  the  design  and  implementation  of  a  new  program  for  the  display 
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and  manipulation  of  waveforms.  The  previous  program  was  originally  designed  for  the 
Tektronix  4014,  a  storage  display  device.  This  program  was  then  converted  to  be  used  on  the 
Megatek  7000,  a  refresh  display  device.  We  are  now  working  on  a  new  program  which  takes 
advantage  of  the  capabilities  of  the  Megatek  terminal,  which  include  the  ability  to  segment  the 
elements  on  the  display  screen  and  manipulate  them  individually. 

The  program  accepts  commands  from  disk  files,  the  user's  keyboard,  and  from  menu  items 
on  the  Megatek  screen  selected  by  using  the  data  tablet,  a  peripheral  device  connected  to  the 
Megatek,  which  tracks  the  location  of  a  pen  on  a  1— ft2  tablet.  Figure  1-7  shows  a  typical  display, 
with  the  menu  boxes  on  the  right,  the  waveforms  in  the  center,  and  boxes  describing  the  wave¬ 
form  on  the  left.  Commands  from  all  sources  are  converted  to  a  common  set  of  tokens,  which 
are  then  executed.  Therefore,  additional  command  sources  such  as  a  knob  box  or  a  voice  ana¬ 
lyzer  can  be  interfaced  to  the  system  in  an  orderly  manner.  Commands  are  entered  in  postfix 
notation,  with  the  arguments  preceding  the  command  name.  This  reduces  the  number  of  key¬ 
strokes  required  of  the  user,  and  eliminates  the  need  for  a  delimiter  at  the  end  of  a  variable  set 
of  arguments  and  the  need  for  a  separator  symbol  between  commands. 
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In  addition  to  the  capabilities  available  with  the  previous  display  program,  the  new  program 
allows  the  waveforms  to  be  altered  dynamically,  without  requiring  the  entire  screen  to  be  re¬ 
plotted.  Using  a  joystick  as  an  interface,  the  new  program  allows  smooth  vertical  and  horizontal 
gain  change,  as  well  as  left  and  right  scrolling  of  waveforms.  The  program  allows  phase  names 
to  be  selected  upon  entry  to  the  program,  and  allows  them  to  be  changed  during  the  execution  of 
the  program.  We  plan  to  add  new  functions  to  the  program,  such  as  waveform  filtering  and 
decimation. 

We  have  developed  a  new  set  of  interface  routines  for  use  with  the  Megatek  7000,  to  take 
the  place  of  those  originally  provided  by  Purdue  University.  The  new  routines  follow  the  gen¬ 
eral  outlines  of  the  ACM  Siggraph  Core  Standard  for  graphics  interfaces,  but  are  optimized  for 
both  the  SAS  application  and  the  Megatek  7000. 

The  interface  routines  recognize  waveforms  as  a  special  type  of  graphics  segment  and  ex¬ 
ploit  a  short  vector  format  supplied  by  the  Megatek  to  effectively  double  the  number  of  data 
samples  which  can  be  displayed  without  screen  flicker.  With  the  routines  supplied  by  Purdue, 
noticeable  flicker  developed  with  three  or  more  waveforms  on  the  screen.  With  the  new  rou¬ 
tines,  six  waveforms  can  be  displayed  without  flicker. 

The  new  routines  also  supply  graphics  capabilities  not  found  in  the  Purdue  routines,  most 
importantly  in  the  area  of  dynamic  picture  modification,  as  in  the  case  of  smooth  scrolling  of 
waveforms.  In  addition,  the  routines  make  more  effective  use  of  the  Megatek  display  memory, 
allowing  larger  amounts  of  the  waveform  database  to  be  loaded.  This  results  in  less-frequent 
requests  for  data  from  the  host  processor  disk. 

P.  T.  Crames 
H.  Lison 


G.  ANALYST  INTERACTION  WITH  DATABASE 

Following  Arrivals  Based  Analyst  Review  (ABAR),  the  database  markerfile  contains  adjust¬ 
ments  to  PPK's  automatic  picks,  additional  picks  of  later  phases,  and  deletions  of  phases  con¬ 
sidered  spurious  by  the  analyst.  It  also  contains  markers  corresponding  to  any  amplitudes  the 
analyst  may  want  to  measure.  These  markers,  together  with  information  from  the  .  gi  file,  are 
input  to  the  program  CPAR.  CPAR  functions  both  as  a  reformatting  program  and  a  signal- 
analysis  program.  It  is  a  reformatting  program  because  it  puts  information  it  finds  in  the 
markerfile  into  INGRES  "arrivals  primitives"  format  and  into  human  readable  format  for  out¬ 
put  to  the  terminal.  It  is  a  signal-analysis  program  because  it  makes  measurements  on  the 
seismic  signal. 

As  a  reformatting  program,  CPAR  uses  the  marker  name  and  sample  number  from  the 
database  markerfile.  The  marker  name  contains  the  phase  name,  INGRES  identification  num¬ 
ber,  and  information  which  identifies  it  as  an  arrival-time  or  amplitude  marker.  The  sample 
number  is  used  to  calculate  the  phase  arrival-time  or  to  locate  the  point  on  the  seismic  trace 
to  be  analyzed. 

The  program  utilizes  marker-name  information  in  the  following  way.  An  "A"  in  the  first 
column  designates  the  marker  as  an  amplitude.  Any  other  letter  designates  it  as  an  arrival 
marker  name.  For  markers  designated  as  corresponding  to  arrivals,  there  are  several  letters 
which  have  special  meanings  when  placed  in  the  first  column.  A  "D"  in  the  first  column  signifies 
an  analyst  decision  during  the  review  session  that  the  arrival  is  either  misidentified  or  false. 
CPAR  creates  a  "remove  arrival"  primitive  for  any  arrival  whose  marker  name  begins  with 
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a  "D."  Other  special  letters  are  "i,"  "e,"  and  *w.*  CPAR  understands  these  as  estimates  of 
arrival-time  reading  accuracy.  The  analyst  calls  the  arrival  *1"  (impulsive),  if  he  estimates 
the  accuracy  of  his  pick  to  he  better  than  0.2  s.  For  an  estimated  accuracy  of  0.2  to  1.0  s,  the 
arrival  is  called  "e"  (emergent).  A  *w*  (weak)  in  the  first  column  designates  an  arrival-time 
accuracy  of  worse  than  1  s. 

CPAR  searches  columns  two  through  seven  of  both  arrival  and  amplitude  marker  names, 
looking  for  the  phase  name  of  the  arrival.  A  *t"  occurring  in  columns  three  or  four  of  an 
arrival-time  marker  means  that  the  associated  event  has  been  identified  as  a  teleseism.  An 
"1"  found  in  these  columns  categorizes  the  event  as  a  local.  CPAR  inspects  the  phase  name  of 
an  amplitude  marker  in  order  to  associate  the  ground  motion  and  period  it  calculates  with  the 
proper  arrival.  The  first  blank  encountered  in  a  column  prior  to  the  seventh  terminates  the 
phase  name. 

Starting  with  the  first  column  following  the  blank  column  in  the  arrival  marker  name,  CPAR 
looks  for  the  INGRES  identification  number,  a  unique  number  of  up  to  5  digits  assigned  to  each 
arrival  in  the  INGRES  database.  If  it  finds  the  INGRES  ID  number,  CPAR  creates  an  "update 
arrival"  primitive  for  those  arrivals  not  previously  eliminated  because  of  a  "D"  in  the  first 
column  of  the  marker  name.  If  no  ID  number  is  found,  the  arrival  is  assumed  to  have  been 
added  during  the  ABAR  session.  In  this  case,  CPAR  places  the  results  of  its  analysis  in  a 
"new  arrival"  primitive. 

Arrival  times  are  found  by  combining  the  seismogram  start  time  and  sample  rate,  supplied 
from  the  .gi  file,  with  the  sample  number  of  the  arrival  found  in  the  markerfile.  Times  are 
written  to  the  terminal  in  the  standard  format:  month/day/year  hour:minute:second. tenths. 

The  format  for  the  INGRES  "arrival  primitive"  is  hard  for  a  human  to  decipher.  Year  and  day 
of  year  are  represented  by  a  5-digit  integer,  where  the  first  2  digits  on  the  left  represent  the 
year  and  the  last  3  represent  the  day  of  the  year.  The  time  of  day  is  in  seconds.  CPAR  also 
finds  the  difference  between  the  first  arrival  and  all  succeeding  phases  on  the  seismogram. 

This  information  is  written  to  the  terminal  for  the  convenience  of  the  analyst. 

CPAR  performs  several  analysis  functions.  It  finds  the  direction  of  first  motion  of  any 
first  arrival  which  the  analyst  has  called  an  "i."  For  any  signal  categorized  as  teleseismic,  it 
automatically  finds  the  greatest  ground  motion  and  associated  period  in  the  first  5  s  of  the  sig¬ 
nal.  Furthermore,  it  will  find  the  ground  motion  of  any  peak  marked  by  the  analyst,  the  only 
condition  being  that  the  amplitude  marker  name  be  associated  with  an  identical  arrival-time 
phase  name.  Ground  motion  is  one-half  peak-to-peak  in  millimicrons  corrected  for  the  fre¬ 
quency  response  of  the  instrument.  The  period  is  two  times  the  absolute  value  of  the  time 
difference  between  the  measured  peak  (or  trough)  and  the  greatest  adjacent  trough  (or  peak). 

M.  W.  Shields 


H.  AUTOMATIC  ASSOCIATION 

The  automatic  association  process  attempts  to  determine  event  hypocenters  from  the  time- 
ordered  arrival  list.  The  latter  is  obtained  by  merging  the  Level-I  parameters  (produced  by 
the  analyst  at  the  station)  with  those  parameters  produced  as  a  result  of  the  automatic  pro¬ 
cessing  as  described  above,  on  the  continuous  and  segmented  waveform  data  tapes.  We  have 
implemented  the  association  procedure  for  2  October  1980  data.  On  that  day  there  were 
1051  short-period  arrivals  in  the  Level-I  data  collected  by  Sweden,  of  which  565  were  primary 
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(i.  e.,  P  or  PKP).  The  automatic  processing  of  the  23  stations  of  waveform  data  comprising  the 
U.S.  contribution  to  the  data -collection  experiment  produced  648  arrivals.  Of  these,  140  were 
identical  (arrival  time  same  to.$  1  s)  to  arrivals  in  the  Level-I  data,  and  369  were  flagged  as 
being  local,  i.e.,  arrivals  identified  by  the  analysts  as  being  local,  or  with  initial  phases  Pn  or 
Pg,  or  (S-P)  times  less  than  60  s. 

The  automatic  association  scheme  currently  implemented  is  not  yet  fully  incorporated  into 
the  INGRES-based  data  structure.  At  present,  association  is  carried  out  by  a  FORTRAN  pro¬ 
gram  which  operates  upon  the  arrival  list  for  a  specified  time  interval,  normally  one  day,  and 
produces  an  event  bulletin  which  is  then  supplied  to  INGRES  for  refinement  using  the  SAS.  On 
the  latter  subsystem,  arrival  and  event  information  were  up'dated  using  the  INGRES  database. 

The  FORTRAN  association  scheme  consists  basically  of  a  series  of  procedures  of  event 
hypothesizing,  testing,  and  merging.  It  does  not  currently  include  the  standard  association 
procedure,  such  as  that  used  by  the  USGS  National  Earthquake  Information  Service,  which  in¬ 
volves  testing  all  possible  combinations  of  arrivals  as  event  hypotheses.  Nevertheless,  it  is 
both  rapid  (<10  min.  VAX  CPU  time  to  process  24  h  of  arrivals)  and,  we  feel,  surprisingly 
successful  in  view  of  its  simplicity.  The  various  steps  of  the  procedure  follow. 

1.  Initial- Event  Hypothesis 

The  arrivals  list  for  2  October  1980  contains  49  event  locations  contributed  by  seismic  array 
sites  NAO,  SUF,  HFS,  GRFA,  EKA,  GBA,  and  YKA,  and  these  form  the  backbone  of  the  initial- 
event  list.  Next,  all  arrivals  flagged  as  local  are  treated  as  event  hypotheses,  the  event  loca¬ 
tion  being  that  of  the  station  and  the  origin  time  being  set  as  that  of  the  arrival  itself.  Finally, 
groups  of  arrivals  from  four  or  more  close  stations  were  used  to  generate  trial  events  by  fitting 
a  plane  wave,  as  if  the  stations  were  a  conventional  seismic  array.  All  hypotheses  were  assigned 
zero  depth  (surface  focus). 

2.  Initial- Event  Testing 

Each  of  the  hypotheses  generated  above  is  now  simply  tested  by  counting  the  number  of 
arrivals  consistent  with  it.  If  a  given  station  has  an  onset  time  within  80  s  of  that  predicted  by 
the  event  hypotheses,  it  is  considered  to  be  possibly  associated  with  the  event.  There  is  no 
guarantee  that  the  number,  distribution,  and  quality  of  the  observations  are  sufficient  to  locate 
the  event.  Hypotheses  with  less  than  four  consistent  arrivals  were  rejected.  In  order  to  some¬ 
what  reduce  the  number  of  redundant  hypotheses,  arrivals  which  were  deemed  consistent  with 
either  array-  or  local- generated  events  were  temporarily  deleted  from  the  list  from  which  the 
hypotheses  generated  by  the  plane  wave  fit  to  groups  of  close  stations  are  formed. 

3.  Initial  Merging 

At  this  stage  we  certainly  have  many  duplicate  hypotheses  of  the  same  events,  and  it  is  wise 
to  attempt  to  reduce  this  redundancy  before  continuing  the  location  procedure.  Event  hypotheses 
were  deemed  to  be  duplicate  if  their  origin  times  differ  by  less  than  100  s,  and  their  locations 
by  less  than  10".  The  event  with  the  largest  number  of  consistent  arrivals  as  determined  by 
step  2  is  retained. 


4.  Second-Event  Testing 

Each  of  the  remaining  hypotheses  is  now  tested  for  validity  using  a  standard  hypocentral 
location  procedure  (Geiger's  method)  with  the  event  constrained  to  a  surface  focus.  If  the  so¬ 
lution  fails  to  converge  after  a  certain  number  of  iterations  (here  set  to  15)  or  if  there  are  less 
than  4  "defining"  arrivals,  the  hypotheses  are  rejected.  "Defining"  arrivals  were  arbitrarily 
set  to  be  those  with  residuals  of  3.5  s  or  loss. 

5.  Second-Event  Merging 

The  remaining  event  hypotheses  are  now  merged  as  in  step  3,  except  that  events  are  now 
considered  duplicate  if  their  origin  times  differ  by  less  than  60  s  and  their  location  by  less 
than  5°. 

6.  Final  Location  Process 

So  far,  events  have  been  constrained  to  surface  focus.  The  remaining  hypotheses  are  now 
retested  with  the  depth  unconstrained.  If,  during  this  location  procedure,  the  final  depth  be¬ 
comes  greater  than  100  km,  the  location  obtained  is  checked  against  a  stored  list  of  seismic 
regions  in  which  only  shallow  events  have  occurred  historically.  If  the  event  lies  in  one  of  these 
regions,  its  depth  is  reconstrained  to  be  surface  focus. 

Analyst  Refinement:—  The  analyst  now  studies  the  event  solutions  for  inconsistencies 
(e.g.,  in  amplitude)  observed  in  the  parameter  data  and  displays  all  the  waveforms  associated 
with  a  particular  event.  In  the  latter  process,  the  analyst  concentrates  on  identifying  depth 
phases  and,  for  close  arrivals,  consistent  (S-P)  time  differences. 

7.  Results  for  2  October  1980  Data 

Of  the  17  events  given  in  the  NEIS  monthly  summary  for  this  date,  only  11  had  4  or  more 
arrivals  in  the  input  arrival  list  supplied  to  the  automatic  association  procedure  used  above, 
and  the  remaining  6  could  therefore  not  possibly  be  located  on  the  basis  of  the  information  given. 

Sixteen  events  were  obtained  by  the  automatic  association  procedure,  and  these  include  all 
but  1  of  the  locatable  11  events  given  in  the  NEIS  monthly  summary.  The  results  are  shown  in 
Table  1-3,  where  events  common  with  the  monthly  summary  are  identified  by  an  asterisk.  Of 
the  six  new  events,  numbers  5,  6,  and  7  are  almost  certainly  real,  but  we  are  less  confident  in 
numbers  9,  13.  and  14,  which  are  all  in  the  Western  U.S.  Subsequent  analyst  study  of  the  latter 
three,  particularly  of  the  waveforms  corresponding  to  the  associated  arrivals,  cast  serious 
doubt  on  the  validity  of  these  events  generated  by  the  automatic  association  process. 

The  final  column  in  Table  1-3  gives  the  sources  of  the  initial  event  hypotheses  leading  to 
the  16  solutions  given.  Ten  events  were  formulated  on  the  basis  of  two  or  more  array  observa¬ 
tions,  5  on  the  basis  of  two  or  more  arrivals  flagged  as  local,  and  one  by  the  plane-wave  fit 
method  to  groups  of  close  stations.  Of  the  original  1073  input  primary  arrivals,  684  remained 
unassociated;  510  of  these  were  flagged  as  being  local  in  character. 

R.  G.  North 
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TABLE 

1-3 

AUTOMATIC  ASSOCIATION  OUTPUT  FOR  2  OCTOBER  1980 

Origin  Time 

Latitude 

Longitude 

d 

No.  of 

No.  of 

Event 

1  (h:m:s) 

<°N) 

(°E) 

(km) 

Arrivals 

Waveformst 

Source 

01:12:39.7 
03:42:45.9 
03:51:11.0 
05:26:04.6 
06:44:42.9 
08:13:09.2 
11:21:12.4 
15:49:34.5 
16:13:50.7 
17:23:12.6 
17:47:36.1 
19:06  :50.6 
20:48:57.2 
21:21:40.0 
22:00:10.7 
23:08:10.2 


51.69 

50.2! 

-35.43 

-19.27 

-22.24 

52.28 

40.68 

-20.92 

40.09 

52.95 

37.56 

20.02 

44.83 

39.53 

43.33 

38.03 


08.43 
30.33 
178.95 
68.56 
70.16 
52.60 
36.68 
-178.54 
-107.74 
18.06 
141.06 
122.34 
-115.90 
-112.72 
135.32 
30.89 


t  Number  of  arrivals  for  which  waveform  data  are  available. 


7  arrays 
3  arrays 
3  arrays 

2  arrays 

3  locals 
2  arrays 
GROUP 
2  arrays 
2  locals 

4  locals 

5  arrays 

5  arrays 
2  locals 

2  locals 

6  arrays 

3  arrays 
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Fig.  1-2.  Schematic  hardware  configuration  of  Prototype  SDC 


Fig.  1-3.  Structure  of  parametric  data.  Fig.  1-4.  Connections  between  two  hosts. 


FUTURE  FUNCTIONS  HANG 
DIRECTLY  UNDER  THIS 
FUNCTION  UMBRELLm 


Fig.  I-  5.  1 -Mbaud  DMR  link 


Fig.  1-6.  EFTP  (Easy  File 
Transfer  Protocol)  receiver 
and  transmitter. 


rr.  SEISMIC  PROCESSING 


A.  GLITCH  NOISE  RECOGNITION  USING  AUGMENTED  TRANSITION 
NETWORK  GRAMMARS 

Nonseismic  signals  such  as  calibration  sequences,  isolated  glitches,  glitch  sequences 
(called  buzzes),  and  boxes  commonly  cause  automatic  signal  detectors  to  produce  false  alarms. 
Figures  II- 1  to  n-3  show  some  examples.  Thus,  part  of  the  design  of  a  good  seismic  detector 
must  be  devoted  to  designing  a  glitch  noise  recognizer.  Running  median  filters1  can  identify 
narrow  isolated  glitches,  but  may  miss  more  complicated  noise  patterns.  Ad  hoc  post-detection 
logic  is  often  added  to  recognize  glitches  by  features  such  as  their  duration  or  one-sidedness. 
Such  tests  can  be  difficult  to  implement  and  may  not  be  completely  effective. 

On  the  other  hand,  a  simple  Augmented  Transition  Network  (ATN)  grammar  has  been  de¬ 
signed  that  effectively  recognizes  various  forms  of  glitch  noise.  The  grammar  can  not  only 
significantly  reduce  the  false-alarm  rate  of  a  signal  detector,  but  can  also  produce  English  de¬ 
scriptions  of  noise  waveforms. 

Since  ATN  grammars  have  been  introduced  in  a  previous  SATS,  they  will  be  described  only 
briefly  here.  A  waveform  grammar  is  a  set  of  rules  that  can  be  used  to  describe  a  waveform. 
For  example,  consider  the  calibration  sequence  shown  in  Fig.  II- i.  A  grammar  rule  that  de¬ 
scribes  a  calibration  sequence  might  be: 

A  positive  pulse  followed  by  60  s  of  noise  followed  by  a  negative  pulse. 

Additional  rules  would  further  describe  the  three  components  of  the  calibration  pulse  (positive 

pulse,  noise,  and  negative  pulse).  The  ATN  formalism  provides  a  powerful  programming  lan- 

3 

guage  in  which  grammar  rules  can  be  expressed  easily. 

A  glitch  recognition  system  for  the  station  ADK  has  been  designed  that  operates  in  five 
steps: 

St:  Detect  possible  seismic  signal 

S2:  Hysteresis  filter  of  detected  waveform 

S3:  Segment  filtered  data  into  words 

S4:  Apply  ATN  grammar  to  words  to  produce  a  description  of  the 
waveform 

S5:  Decide  if  detection  is  due  to  glitch  noise. 

Si:—  The  seismic  detector  used  has  been  described  in  a  previous  SATS.2  This  detector  does 
well  when  compared  with  a  human  analyst,  in  fact,  it  detects  weak  signals  more  reliably  than 
an  analyst  does.  Unfortunately,  when  used  on  noisy  data  from  ADK,  many  false  alarms  are 
produced.  The  detector  identifies  a  time  interval  beginning  60  s  before  the  detection  and  ending 
at  a  time  that  depends  on  the  characteristics  of  the  signal. 

S2:—  After  the  detection  stage,  the  raw  data  are  hysteresis  filtered.  The  output  of  a  hyster¬ 
esis  filter  h(i)  at  time  i  is. 

h(0)  =  x(0)  ; 

h(i)  =  x(i)  ,  |x(i)  —  x(i  —  1)  |  >  c  ; 

=  h(i  —  1)  ,  otherwise  ; 
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where  x(i)  is  the  input  time  series,  and  c  is  a  threshold  value.  Thus,  a  hysteresis  filter  re¬ 
turns  the  value  of  the  input  point  if  it  differs  significantly  from  the  last  point;  otherwise,  it  out¬ 
puts  the  previous  output.  The  dotted  line  in  Fig.  II-4  shows  the  effect  of  hysteresis  filtering  a 
calibration  pulse.  The  effect  of  hysteresis  filtering  is  to  pass  large  variations  unchanged  while 
reducing  low-amplitude  background  noise  to  flat  segments  separated  by  short  jumps.  This  com¬ 
presses  the  data  and  simplifies  the  grammar. 

S3;—  After  hysteresis  filtering,  the  data  are  divided  into  segments  having  either  positive, 
negative,  or  zero  slope  by  a  simple  ATN  grammar.  These  segments,  represented  by  their 
endpoints,  are  the  words  used  as  input  to  the  glitch  grammar. 

S4:—  The  glitch  grammar  was  designed  on  18  data  samples  from  ADK  recorded  on  4  Octo¬ 
ber  1980.  Only  one  of  these  signals  was  a  true  seismic  signal.  Since  hysteresis  filtering  tends 
to  destroy  weak  seismic  signals,  no  attempt  was  made  to  identify  true  seismic  signals.  The 
only  requirement  was  that  the  glitch  grammar  passed  true  signals  without  declaring  any  part  of 
them  to  be  glitch  noise. 
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The  glitch  grammar  produces  a  structural  description  of  the  time  series.  This  description 
can  be  translated  by  computer  into  an  English  description  that  is  similar  to  a  seismologist's. 
Such  descriptions  are  used  to  caption  Figs.  II- 1,  II-2,  and  n-3. 

S  5.-—  Based  on  the  description  produced  by  the  grammar,  each  detection  was  assigned  to 
one  of  two  classes,  S  or  G.  The  signal  was  assigned  to  class  G  (glitch  noise)  if  a  box,  glitch, 
buzz,  or  calibration  pulse  occurred  within  2  s  of  the  detector  onset  time;  otherwise,  the  detec¬ 
tion  was  assigned  to  class  S  (possible  seismic  signal). 

This  detector  was  applied  to  24  h  of  data  recorded  at  ADK  on  2  October  1980;  39  detections 
occurred.  Table  II- 1  shows  a  confusion  matrix  of  the  results.  Over  7  5  percent  of  the  signals 
are  correctly  classified.  Since  no  class  S  signal  was  misclassified,  the  detector  is  a  reliable 
but  conservative  glitch  recognizer. 


TABLE  11-1 

CLASSIFICATION  PRODUCED 

BY  ATN  CLASSIFIER 

Assigned  Class 

Actual  Class 

Seismic 

Glitch 

Seismic 

11 

0 

Glitch 

6 

22 

The  six  class  G  signals  that  were  missed  were  small  glitches  isolated  in  background  noise 
that  did  not  occur  in  the  training  set.  Thus,  the  grammar  is  completely  effective  at  recogniz¬ 
ing  the  class  of  signals  it  was  trained  on  and  can  be  extended  easily  to  handle  this  new  type  of 
glitch.3 


Although  the  class  of  glitches  that  the  grammar  can  recognize  is  quite  large,  no  attempt  was 
made  to  make  it  general  enough  to  work  at  all  stations.  The  reason  for  this  is  that  it  is  better 
to  write  grammars  for  individual  stations  than  to  try  to  find  a  general  theory  of  glitches.  A 
grammar  should  be  thought  of  as  a  local  expert  that  uses  detailed  knowledge  to  recognize  glitches 


at  one  station. 


K.  R.  Anderson 


B.  DYNAMIC  WAVEFORM  MATCHING 

Although  the  glitch  recognition  grammar  presented  above  is  simple,  consisting  of  only 

41  states,  other  seismic  applications  can  require  more  complicated  grammars.  For  example, 

4 

I.iu  and  Fu  developed  a  grammar  for  distinguishing  between  earthquake  and  nuclear-explosion 
signals  that  required  nearly  800  states.  Thus,  techniques  that  help  automate  the  grammar  de¬ 
sign  process  are  valuable.  To  this  end,  the  technique  of  Dynamic  Time  Warping  (DTW)  is  being 
investigated. 

DTW  is  a  technique  used  in  speech  processing  to  match  two  signals  whose  features  need  not 
be  perfectly  time  aligned.  This  technique  is  used,  for  example,  to  match  a  word  utterance  to 
a  set  of  word  templates.  DTW  does  a  nonlinear  time  alignment  between  two  patterns  based  on 
a  least-error  criterion.  In  this  way,  it  compensates  for  temporal  variations  between  the  two 
patterns  due  to  differences  in  the  way  a  speaker  says  the  same  utterance  on  different  occasions. 

We  plan  to  exploit  the  similarity  between  DTW  and  error-correcting  parsing.^  Both  DTW 
and  parsers  match  a  description  of  a  signal  against  an  input  signal.  In  parsing,  the  description 
is  in  the  form  of  grammar  rules.  In  DTW,  the  description  is  in  the  form  of  a  matching  function 
and  a  template  (another  signal). 

Although  grammatical  techniques  are  more  general  than  DTW,  they  are  harder  to  develop. 
Thus,  we  will  use  DTW  as  a  tool  to  understand  seismic  waveform  recognition  problems.  It  is 
not  clear  how  far  DTW  will  take  us  before  we  must  use  more  grammatical  methods,  but  there 
is  a  range  of  methods  between  these  two  extremes.  For  example,  a  template  matcher  that  uses 
local  constraints  on  waveform  shape  behaves  like  a  simple  parser. 

DTW  is  promising  in  several  recognition  tasks  involving  the  envelopes  of  local  seismic 
signals : 

(1)  Template  matching  of  signals  with  known  characteristics  against  new 
incoming  signals.  For  example,  template  matching  could  be  used  to 
identify  portions  of  the  signal  such  as  the  P- arrival,  the  S-ar rival,  and 
the  coda.  The  interval  between  the  P-  and  S-arrivals  can  be  used  to 
estimate  the  distance  of  the  event,  and  information  from  the  coda  can 
be  used  to  estimate  spectral  parameters  of  the  earthquake  as  well  as 
regional  variations  in  attenuation.  Some  seismic  stations  receive 
most  of  their  local  earthquake  activity  from  a  small  number  of  local 
source  regions,  such  as  quarrys.  Template  matching  could  be  used 

to  identify  signals  from  these  source  regions. 

(2)  Cluster  analysis  of  seismic  envelopes  to  aid  in  development  of  tem¬ 
plate  matching  techniques  and  as  an  initial  step  in  grammatical 

;  ,  4.7 
inference. 


After  a  brief  introduction  to  DTW,  progress  in  each  of  these  areas  is  summarized. 


1.  Dynamic  Time  Warping  (DTW) 

DTW  can  be  described  as  follows:  Given  a  reference  pattern  R(n),  1  <  N,  and  a  test 
pattern  T(m),  4m^M,  find  an  optimal  mapping  from  each  time  axis  n  and  m  onto  a  common 
time  axis  k: 

n  =  i(k)  ,  k  =  1,  2.....K 

m  =  j(k)  ,  k  =  1,2,...,K  . 

In  general,  R(n)  and  T(n)  can  be  multidimensional  feature  vectors,  but  in  this  study  they  are 
scalar  quantities  representing  the  signal  envelope.  The  optimal  mapping  is  one  that  minimizes 
a  total  distance  function  D  of  the  form 

K 

D  =  min  4  2  d{Rfi(k)),T[j(k)}>  w(k)  (II- 1) 

[i(k),j(k))  w  k=1 

where  d[R(i),T(j)l  is  the  local  distance  between  the  point  i  of  the  reference  and  point  j  of  the 
template,  w(k)  is  a  local  weighting  function,  and  W  is  a  corresponding  normalization.  Several 
forms  of  Eq.  (II— 1 )  have  been  proposed  and  their  performance  in  word- recognition  problems  has 
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been  studied  by  Myers. 

It  is  convenient  to  think  of  minimizing  Eq.  (II- 1)  in  terms  of  finding  an  optimum  path 
(i(k),j(k))  in  (n,  m)  space  as  shown  in  Fig.  n-5.  To  solve  Eq.  ( II- 1)  as  a  path  finding  problem, 
one  must  specify  several  things: 

(a)  Endpoint  constraints  —  the  way  that  the  path  begins  and  ends 

(b)  Local  continuity  constraints  —  constraints  on  motion  from  one  point  to 
another 

(c)  Global  path  constraints  —  limitations  on  where  the  path  can  wander  in 
(n,  m)  space 

(d)  Range  constraints  —  the  maximum  separation  between  i(k)  and  j(k)  in  a 
valid  match 

(e)  Axis  orientation  —  the  effects  of  interchanging  the  roles  between  the 
test  and  reference 

(f)  Distance  measures  —  both  the  local  measure  of  similarity  and  the  over¬ 
all  distance  function  used  to  determine  the  optimal  path. 

Myers  discusses  each  of  these  factors  and  evaluates  the  performance  of  several  combina¬ 
tions  of  them.  We  have  investigated  several  of  these  alternatives  for  seismic  envelope  matching 
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and  have  found  that  his  type  II-d  matching  function  works  reasonably  well. 

The  solution  to  Eq.  (II- 1)  can  be  written  as 

D(N,  M) 

D  =  -  (II-2) 

where  D_  is  called  an  accumulated  distance  function,  which  for  the  case  considered  here  is  re- 
a  ' 

cursively  defined  as 
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Da(n,  m)  =  min 


[  Da(n~  1,  m-  i)  +  2d(n,  m)  "| 


Da(n-l,m-2)  +  3d(n,  m) 


Da(n-2,  m-1)  +  3d(n,  m) 


with  the  initial  condition 


Da(l,i)  =  d(l,l) 


This  type  of  recursive  solution  to  Eq.  ( II- 1)  is  known  as  Dynamic  Programming. 


2.  Waveform  Matching 

Figure  11-6  shows  an  example  of  the  DTW  of  the  envelopes  of  two  local  seismograms.  The 
envelope  is  constructed  by  computing  nonoverlapping  3-s  block  averages  of  the  rectified  wave¬ 
forms.  This  averaging  smooths  out  high-frequency  fluctuations  in  the  envelope  but  provides 
enough  detail  for  this  application.  Before  matching,  the  envelopes  are  normalized  by  their  max¬ 
imum  amplitude.  The  local  matching  function  used  is: 

d[R(i),T(j)]  =  |R(i)  -T(j)|  . 

In  other  words,  the  waveforms  are  aligned  based  on  the  area  under  their  normalized  envelopes. 

Dotted  lines  in  the  figure  are  drawn  between  a  point  in  each  waveform  that  corresponds  in 
the  match.  Note  that  the  onsets  of  the  P-  and  S-arrivals  have  been  matched. 

Figure  II- 7  shows  the  match  between  two  waveforms  that  are  not  so  similar.  The  lower 
envelope  has  a  much  broader  P-arrival  than  the  top  and  no  correspondence  is  found.  Also,  the 
correspondence  between  the  S-arrivals  is  not  very  good  due  to  the  difference  in  their  shape. 

Thus,  time  warping  might  be  used  to  identify  P-  and  S-arrivals  by  matching  an  unknown  envelope 
against  a  set  of  templates  and  measuring  the  P-S  interval  from  the  best  match,  but  care  must 
be  taken  in  choosing  the  set  of  templates. 


3.  Waveform  Clustering 

Figure  II-8  shows  how  DTW  can  be  used  to  cluster  waveforms.  The  right  part  of  the  diagram 
shows  the  envelopes  of  eight  local  seismograms  from  station  BDW.  In  the  left  part  of  the  figure, 
a  dendrogram  is  used  to  show  the  hierarchical  clustering  of  the  waveforms.*^  A  dendrogram  is 
a  tree  that  shows  the  hierarchical  clustering  of  a  set  of  samples  (in  this  case,  waveforms).  The 
leaf  nodes  of  the  tree  are  the  waveforms  themselves.  Two  nodes  are  merged  into  one  at  a  level 
indicating  their  dissimilarity.  For  example,  the  dendrogram  shows  that: 

(a)  Waveforms  1  and  2  are  most  similar,  followed  by  5  and  6. 

(b)  Waveforms  i  through  4  are  relatively  similar  to  each  other,  as  are 
waveforms  5  through  7. 

(c)  Waveform  8  is  most  unlike  the  other  waveforms. 

This  clustering  is  similar  to  how  a  seismologist  would  cluster  these  waveforms. 

The  next  step  in  the  research  is  to  cluster  a  large  number  of  waveforms  from  several  sta¬ 
tions.  These  clusters  will  be  used  as  templates  in  the  design  of  a  template  matcher  and  as  sen¬ 
tences  in  the  design  of  waveform  grammars. 

K.  R.  Anderson 
J.  E.  Gaby 
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C.  WAVEFORM  SEGMENTATION  USING  AFFINITY 

The  first  step  in  syntactic  waveform  analysis  is  to  divide  the  waveform  into  segments  of 
morphological  significance  called  words.  This  sequence  of  words  should  be  a  simplified  repre¬ 
sentation  of  the  waveform,  but  must  have  enough  information  to  extract  features  of  interest. 

The  "affinity"  technique*1  that  has  been  used  for  region  growing  in  scene  analysis  is  being  used 
to  produce  word  sequences  from  the  envelopes  of  local  and  regional  seismograms. 

Affinity  is  a  general  technique  that  organizes  data  into  a  binary  tree  structure.  It  groups 
together  regions  based  on  their  relative  similarity.  Two  similar  regions  are  merged  together 
to  form  one  region  at  a  higher  node  in  the  tree.  The  new  node  is  an  abstracted  representation 
of  the  two  regions  below  it. 

There  are  several  reasons  why  affinity  is  an  attractive  approach.  First,  the  similarity 
criterion  is  relative,  no  thresholds  need  be  set;  the  waveform  essentially  segments  itself.  Sec¬ 
ond,  it  preserves  the  information  content  of  the  signal.  To  acquire  more  detailed  information 
about  the  signal,  one  need  only  go  deeper  down  the  tree.  Although  the  leaves  of  the  tree  could 
be  the  individual  data  points  themselves,  it  is  more  efficient  to  cut  the  tree  off  at  a  level  that 
provides  adequate  resolution.  Third,  the  tree  is  a  flexible  data  structure  that  allows  the  wave¬ 
form  to  be  divided  into  word  sequences  in  several  ways.  Fourth,  a  grammar  may  parse  a  par¬ 
ticular  sequence  of  words  derived  from  a  certain  level  in  the  tree.  Information  at  lower  levels 
in  the  tree  can  then  be  used  to  refine  or  verify  the  parse  or  extract  useful  features. 

1.  Introduction  to  Affinity 

The  affinity  algorithm  compares  a  segment  of  a  waveform  with  its  adjacent  neighbors.  Seg¬ 
ments,  whose  descriptions  closely  match,  are  then  merged  together  to  form  a  new  segment  that 
is  representative  of  the  two  similar  segments.  The  ith  segment  is  described  by  the  feature  vec¬ 
tor  x(i),  chosen  for  the  particular  application. 

Affinity  compares  one  segment  x(i)  to  its  adjoining  neighbors  x(i  —  1)  and  x(i  +  1).  A  dis¬ 
tance  function  F  determines  if  x(i)  is  more  similar  to  x(i  —  1)  or  x(i  +  1).  If  |  F[  x(  i) ,  x(  i  —  1)]|  < 

|  F[x(i),  x(i  +  1)|  |,  segment  i  is  linked  to  segment  i  -  1;  otherwise,  it  is  linked  to  i  +  1.  This 
is  done  for  all  i.  Two  segments  (a  and  b)  are  merged  if  and  only  if  segment  a  is  linked  to  b 
and  segment  b  is  linked  to  a.  A  merge  function  M[a,  b)  takes  two  segments  and  merges  them 
into  one  segment. 

The  tree  structure  is  developed  from  a  signal  by  successive  applications  of  linking  segments 
of  a  waveform  and  merging  linked  segments.  Segments  are  first  linked  to  their  appropriate  neigh¬ 
bor.  Then  these  linked  segments  are  merged.  This  forms  a  new,  shorter  sequence  of  segments 
that  represents  the  same  waveform.  This  process  is  repeated  until  the  entire  waveform  has  been 
merged  into  one  segment  whose  set  of  features  represents  the  waveform  as  a  whole. 

The  one  final  segment  represents  the  top  node  of  the  tree.  It  branches  down  to  two  segments 
that  were  merged  to  form  the  top  node.  Likewise,  the  two  branches  are  nodes  that  branch  again, 
in  a  binary  fashion.  Each  node  is  a  generalized  representation  of  the  node  below.  As  you  go 
down  the  ti  ee  structure,  there  is  more  detailed  information  about  the  sequence  of  segments  that 
represents  the  signal. 

2.  Application  of  Affinity  to  Seismic  Waveforms 

To  apply  affinity  one  must  determine  the  initial  length  of  a  segment,  the  feature  vector  that 
describes  each  segment,  and  the  functions  F  and  M.  This  decision  is  based  on  the  nature  of  the 
waveform,  the  nature  of  the  features  to  be  ultimately  extracted,  and  the  resolution  required. 
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The  goal  of  this  research  is  to  automatically  extract  information  from  regional  seismograms, 
such  as  the  S-P  interval.  Since  we  are  initially  interested  in  S-P  intervals  larger  than  5  s, 
a  1-s  initial  segment  size  (20  samples)  provides  adequate  time  resolution. 

In  this  study  three  features  were  analyzed:  the  average  magnitude  (absolute  amplitude),  the 
average  zero-crossing  rate,  and  the  slope  between  two  segments.  These  features  could  be  com¬ 
bined  in  different  ways,  but  they  were  analyzed  separately.  Each  of  these  features  uses  differ¬ 
ent  linking  and  merging  functions. 

Affinity  based  on  average  magnitude  is  used  because  the  envelope  contains  visual  cues  for 
the  P  and  S  onset  times.  Each  segment  is  represented  by  three  parameters  —  its  length,  its 
average  magnitude,  and  the  time  of  the  middle  of  the  segment,  L(i),  A(i),  and  T(i),  respectively. 
The  distance  function  F  is  the  difference  in  the  magnitude  values  between  a  segment  and  its  two 
neighbors  weighted  by  the  segment  lengths 

-,u.,  A(a)L(a)  -  A(b)L(b) 

F[x(a),x(b)| - L(a]  -TLfE) -  * 

Merging  two  segments  is  accomplished  by  producing  a  new  length,  magnitude,  and  time  to 
describe  the  new  segment.  The  new  length  is  the  sum  of  the  two  lengths  of  the  segments  being 
merged.  The  new  magnitude  is  the  weighted  sum  of  the  average  magnitudes  of  the  two  segments 
being  merged.  The  new  time  is  simply  the  center  of  the  new  segment: 

A(new)  -  A(a)L(a)  +  A(b)L(b) 

A(new)  -  L(a)  +  L(b) 

L(new)  =  L(a)  +  L(b) 

T(new)  -  T(a)  +  . 

Affinity  based  on  the  average  number  of  zero  crossings  is  similar  to  the  average  magnitude 
version.  The  only  difference  is  that  A(i)  is  replaced  by  Z(i) ,  the  average  zero-crossing  rate  in 
that  interval.  The  zero-crossing  rate  is  a  rough  measure  of  the  frequency  content  of  the  signal. 

It  can  be  useful  in  identifying  weak  arrivals  buried  in  noise,  because  a  sudden  change  in  the 
zero-crossing  rate  may  be  the  only  evidence  for  a  weak  arrival. 

Affinity  based  on  slope  requires  a  different  set  of  parameters  to  describe  a  segment.  Here, 
the  slope  between  the  average  magnitude  of  two  segments,  of  the  initial  i-s  segmentation,  is 
used  as  the  feature.  Each  slope  segment  is  represented  by  four  numbers  [xl(i),  yl(i),  x2(i),  and 
y2(i))  where  the  x's  are  the  time  index  and  the  y's  are  the  average  amplitude.  The  information 
in  each  segment  is  duplicated  in  the  following  segment.  That  is  to  say,  x2(i)  =  xt(i  +  1)  and 
y2(i)  =  yi(i  +  i).  Because  of  this  overlap,  the  information  content  of  this  initial  sequence  is 
similar  to  that  of  the  average  magnitude  affinity. 

Slope  affinity  compares  the  slopes  of  adjacent  segments  and  links  to  the  segment  with  the 
minimum  change  of  slope.  The  segment  is  replaced  by  a  new  segment  with  the  new  start  point 
the  same  as  the  x-y  coordinates  of  the  start  point  of  the  first  segment,  and  the  endpoint  the 
same  x-y  coordinates  of  the  endpoint  of  the  second  segment.  The  distance  function  Ffa,  b)  is 

F[x(a),  x(b)]  =  yl(a)  +  y2(a)  ♦  y2(b)  »ggj  I  *f|fj 
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and  the  parameters  of  a  new  merged  segment  are 

xl(new)  =  xl(a) 
yl(new)  =  yl(a) 
x2(new)  =  x2(b) 
y2(  new)  =  y2(b) 

3.  Results 

Of  the  three  implementations  studied  (average  magnitude,  slope,  and  zero-crossing  rate), 
the  first  two  yielded  good  results.  The  zero-crossing  density  affinity  creates  word  sequences 
that  can  be  either  misleading  or  show  no  significant  information.  It  provides  some  information 
when  the  P-arrival  is  buried  in  the  noise  level,  but  it  is  not  always  reliable.  A  nonlinear 
F  function  that  uses  magnitude  when  the  signal  was  strong  and  frequency  when  the  signal  was 
weak  (much  as  a  seismologist  does),  might  provide  better  segmentations. 

Affinity  based  on  average  magnitude  produces  more  meaningful  segmentations.  Figure  II-9 
shows  the  tree  structure  that  resulted  from  the  application  of  affinity  on  one  waveform.  In  this 
figure,  and  in  all  others,  the  initial  segmentation  is  based  on  3  s  of  the  waveform  and  not  1  s, 
as  used  in  the  study.  This  makes  the  results  easier  to  see,  as  the  1-s  data  are  highly  detailed 
and  confusing.  Regions  with  a  high  degree  of  similarity  are  correctly  merged  to  form  larger 
segments.  Figure  Il-lO(a-c)  shows  three  different  sequen  ;es  of  words  derived  from  different 
levels  in  the  tree.  The  dotted  lines  are  the  original  sequences  given  to  the  affinity  routine. 

The  bold  lines  are  the  sequences  of  segments  at  different  levels  in  the  tree,  the  zeroth  level 
being  the  top  node  of  the  tree.  The  two  major  peaks  are  the  P-  and  S-waves.  Word  sequences 
produced  from  levels  5  and  7  show  the  arrival  onsets  clearly  without  showing  too  much  detail 
in  the  coda.  Thus,  they  provide  reasonable  segmentations. 

Figure  II- 11  shows  the  tree  produced  by  the  slope  affinity  algorithm.  The  trees  shown  in 
this  figure  and  in  Fig.  11-9  are  generally  similar.  Figure  II-12(a-c)  shows  different  potential 
word  sequences  derived  from  different  depths  of  the  tree.  In  this  segmentation,  peaks  stand 
out  immediately.  Both  magnitude  affinity  and  slope  affinity  reduce  the  waveform's  complexity 
and  maintain  information  necessary  for  the  recognition  process. 

Figure  II-13(a-d)  shows  a  second  method  of  producing  word  sequences  from  a  tree.  Here 
the  word  sequences  are  produced  using  a  maximum  length  word  rather  than  a  fixed  depth  in  the 
tree.  In  (a)  and  (b),  all  word  segments  are  either  less  than  18  s  in  duration  or  are  a  leaf  node 
of  the  tree.  In  (c)  and  (d),  the  length  is  12  s  or  less.  This  shows  how  flexible  the  tree  struc¬ 
ture  is  in  providing  different  word  segmentations. 

The  affinity  algorithm  has  some  drawbacks.  The  merging  does  not  consider  the  relative 
nature  of  merging  things  very  similar  and  those  somewhat  similar  at  the  same  time.  In  other 
words,  the  branches  of  the  tree  all  have  the  same  relative  length.  This  problem  can  be  corrected 
by  not  only  specifying  nodes  and  branches,  but  also  the  relative  lengths  of  the  branches.  Fig¬ 
ure  II-14(a-c)  illustrates  a  second  problem.  The  P-wave  onset  is  just  visible  as  a  step  before 
the  predominate  peak  S-wave.  The  representation  of  the  noise  before  the  P-arrival  is  more 
detailed  at  a  given  level  in  the  tree  than  the  noise  after  the  signal.  This  is  because  the  length 
of  the  noise  before  the  signal  is  less  than  the  noise  after,  and  thus  requires  a  shallower  subtree 
to  represent  it. 


4.  Conclusion 


Studying  the  morphological  structure  of  seismic  signals  is  greatly  aided  by  the  affinity  algo¬ 


rithm.  Affinity  based  on  average  magnitude  or  slope  yields  promising  results.  In  both  cases, 
the  waveform  is  broken  into  sections  that  reveal  the  information  required  to  determine  the  P-S 
time  interval  while  removing  unnecessary  detail. 

Affinity  is  also  useful  for  analyzing  the  waveform's  structure.  Structuring  the  data  in  a 
tree  format  allows  the  waveform  to  be  looked  at  from  different  points  of  view.  The  tree  itself 
could  be  the  structure  given  to  the  grammar  to  analyze.  The  grammar  could  be  smart  enough 
to  use  the  tree,  with  all  the  diversity,  to  break  the  waveform  up  as  it  needs  to  recognize  what 


it  wants. 


J.  E.  Gaby 

K.  R.  Anderson 


D.  AUTOMATIC  ASSOCIATION  USING  EXPERT  SYSTEMS  TECHNIQUES 

An  automatic  association  (AA)  program  aids  seismologists  in  associating  arrivals  with 

12  13 

earthquakes.  Although  several  AA  programs  have  been  developed,  ’  there  are  still  many 
issues  that  need  to  be  worked  on,  among  which  are. 

(1)  AA  programs  can  correctly  associate  about  70  to  80  percent  of  the 
information  available  to  them.  Can  this  success  rate  be  improved? 

(2)  False  events  can  be  created  by  splitting  one  event  into  several  smaller 
ones,  or  merging  arrivals  from  two  events  into  one.  How  can  this  be 
avoided? 


(3)  Many  events  created  by  these  algorithms  are  based  on  only  a  few  ar¬ 
rivals.  Is  there  any  way  to  estimate  the  probability  that  these  events 
really  occurred? 

(4)  What  is  the  best  way  for  an  analyst  to  interact  with  the  A  A  program? 

AX,  an  Expert  System  to  aid  seismologists  in  the  association  problem,  is  under  development. 
It  is  designed  as  a  research  environment  that  can  be  used  to  resolve  these  issues.  AX  differs 
from  other  association  programs  because  it  is  based  on  a  rule  driven  system  capable  of  using 
knowledge  about  the  association  problem  in  much  the  same  way  a  seismologist  would.  Since 
Expert  Systems  ideas  are  new  to  seismology,  this  report  serves  as  an  introduction  to  Expert 
Systems  and  as  a  progress  report  on  AX. 

Expert  Systems  is  a  branch  of  Computer  Science  research  that  studies  ways  to  represent 

human  expertise  in  computer  systems.  In  the  past  few  years,  several  Expert  Systems  have  been 
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developed  as  can  be  seen  in  Table  U-2.  Although  Expert  Systems  usually  do  better  than  an 

average  person,  they  may  not  work  as  well  as  a  true  expert.  PROSPECTOR  works  well  enough 
that  it  can  be  used  as  a  training  tool  for  new  human  mineral  prospecting  experts,  and  even  teaches 
established  experts  something  new. 

Two  basic  problems  in  building  an  Expert  System  are  how  to  organize  and  represent  knowl¬ 
edge  and  how  to  discover  the  rules  required  to  manipulate  that  knowledge  in  the  particular  cog¬ 
nitive  task.  The  appropriate  knowledge  is  usually  uncovered  by  protocol  analysis.  This  involves 
observing  a  human  expert,  such  as  a  seismic  analyst,  doing  a  task  and  having  him  explain  as 
completely  as  possible  his  decisions  and  the  reasoning  behind  them. 


TABLE  11-2 

EXPERT  SYSTEMS  STUDIES 


Program 

Domain 

Reference 

MYCIN 

Medical  diagnosis 

14 

PROSPECTOR 

Mineral  prospecting 

15 

DIPMETER  ADVISOR 

Oil-well  logging 

16 

DENDRAL 

Chemical  analysis 

17 

There  has  been  one  protocol  session  with  a  seismic  analyst  using  an  interactive  association  1 

system  at  SDAC.  A  detailed  log  of  his  interaction  with  the  system  was  kept  as  he  associated  ] 

arrivals  for  several  events.  To  give  some  idea  of  the  work  involved  in  associating  an  event,  the  1 

analyst  went  through  over  50  different  steps  during  the  association  of  a  medium-size  event.  \ 

Additional  protocol  sessions  will  be  necessary  during  the  development  of  AX.  One  advantage  of 
the  Expert  Systems  approach  is  that  AX  can  be  developed  incrementally  based  on  these  protocol 
sessions. 

Once  the  knowledge  used  by  an  analyst  is  available,  it  must  be  represented  in  the  computer. 

Of  course,  any  program  that  solves  a  problem  can  be  said  to  embody  knowledge  of  a  domain. 

The  distinguishing  characteristic  of  an  Expert  System  is  that  much  of  its  knowledge  about  the 
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task  domain  is  stored  explicitly  in  a  database.  AX  uses  three  Expert  Systems  programming 
techniques  to  represent  knowledge: 

(1)  FRL,  a  knowledge  representation  language 

(2)  Sentinels 

( 3)  Rules. 

Each  of  these  techniques  will  be  described  in  turn, 

1.  Frame  Representation  Language  (FRL) 

The  FRL  is  a  package  of  LISP  functions  for  manipulating  data  structures  known  as 
frames.1^’20  Originally  motivated  by  Minsky,^1  frames  have  several  features  that  make  them 
convenient  for  knowledge  representation.  The  objects  used  by  AX  are  represented  as  frames 
in  the  frame  database.  This  includes  seismological  objects  such  as  events  and  arrivals,  and 
also  the  sentinels  and  rules  that  manipulate  them. 

A  frame  is  a  5-level  tree  structure  with  the  following  names  assigned  to  each  level: 

FRAME 

SLOT 

FACET 

DATUM 

COMMENT 

Each  frame  can  contain  several  slots,  each  slot  can  contain  several  facets,  and  so  on  for  each 
level. 

Frames  are  used  to  represent  objects,  slots  are  used  to  represent  the  attributes  of  objects, 
and  the  datum  on  the  value  facet  of  a  slot  represents  the  value  of  the  attribute.  For  example,  an 
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arrival  called  ARRIVAL- IS  that  is  a  PKP  arrival  from  station  ZOBO  occurring  at  time  1234.56 
and  belonging  to  a  related  arrival  group  called  RAG- 12  would  be  represented  as: 

FRAME 
arrival  -  i  6 


SLOT 

FACET 

DATUM 

AKO 

value 

arrival 

time 

value 

12  34.56 

phase 

value 

pkp 

station 

value 

zobo 

rag 

value 

rag- 12 

The  value  of  a  slot  can  be  another  frame,  such  as  ARRIVAL,  PKP,  ZOBO,  and  RAG-12  here. 
The  AKO  (A  Kind  Of)  slot  shows  that  ARRIVAL- 15  is  a  kind  of  ARRIVAL.  The  ARRIVAL  frame 
looks  something  like  this: 

FRAME 

arrival 


SLOT 

FACET 

DATUM 

AKO 

value 

thing 

time 

phase 

required 

(valid -phase  ■ 

value) 

station 

amplitude 

required 

(valid -station 

:  value) 

rag 

slowness 

azimuth 

to-db 

value 

(to-db-arrival 

:frame) 

Most  slots  have  their  obvious  seismological  meaning.  The  TO-DB  slot  describes  how  infor¬ 
mation  in  an  arrival  frame  is  to  be  communicated  to  the  parameter  database.  The  AKO  slot 
allows  frames  to  be  organized  into  an  inheritance  hierarchy.  Properties  of  a  frame  are  inherited 
from  frames  higher  in  the  inheritance  tree.  The  THING  frame  is  at  the  top  of  the  inheritance 
hierarchy.  Thus,  an  ARRIVAL  is  a  kind  of  THING,  and  ARRIVAL- 1  5  is  a  kind  of  ARRIVAL. 

In  addition  to  the  value  facet,  FRL  has  several  facets  that  work  in  combination  with  it 

DEFAULT  A  default  value,  to  be  used  if  no  value  is  available. 

IF-ADDED  A  procedure  executed  if  a  value  is  added  to  a  slot. 

IF'- REMOVED  A  procedure  executed  if  a  value  is  removed  from  a  slot. 

IF -NEEDED  A  procedure  executed  if  a  value  is  needed  but  unavailable. 

REQUIRED  A  check  on  the  validity  of  a  value  added  to  a  slot. 

Such  procedures  work  with  the  inheritance  hierarchy  to  keep  the  frame  database  consistent, 
and  make  deductions.  For  example,  when  ARRIVAL- 1 5  was  created,  checks,  inherited  from 
the  ARRIVAL  frame,  were  made  to  insure  that  the  phase  and  station  name  were  valid. 

The  frames  representing  seismic  parameters  have  an  inheritance  hierarchy  as  shown  in 
Fig.  11-15.  The  EVENT,  RAG,  ARRIVAL,  and  STATION  frames  have  their  obvious  seismic  in¬ 
terpretation.  H- ARRIVALS  are  hypothetical  arrivals  created  by  AX  in  an  attempt  to  correct 
common  arrival  errors.  They  will  be  described  further  below. 

Seismic  waves  are  divided  into  their  own  inheritance  hierarchy  below  the  WAVE  frame  as 
showr  in  Fig.  11-16. 
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2.  Sentinels 

22 

A  sentinel  is  a  procedure  that  waits  for  a  triggering  condition  to  occur  in  the  database, 
and  then  does  some  action.  Sentinels  may  vanish  after  they  complete  their  action,  or  remain 
indefinitely.  For  example,  the  following  sentinel  propagates  information  about  the  distance 
range  of  an  arrival  to  its  related  arrival  group: 

WHENEVER: 

An  arrival  has  a  range  estimate  and  is  assigned  to  a  rag. 

ACTION: 

Assign  its  range  estimate  to  the  rag. 

Sentinels  can  create  other  sentinels  or  groups  of  sentinels.  For  example: 

WHENEVER: 

A  rag  comes  from  a  station  that  has  neighboring  stations. 

ACTION: 

Create  an  event  hypothesis  and  create  sentinels  that 
add  incoming  rags  from  neighboring  stations  to  the  event. 

After  all  arrivals  that  could  be  associated  with  the  event 
have  been  read  in,  remove  all  the  sentinels. 

Thus,  with  sentinels  one  can  create  expectations  of  what  should  happen  to  the  database  in  the 
future,  based  on  its  current  state  of  knowledge. 

3.  Rules 

Typically,  programs  are  specified  in  terms  of  if  —  then  —  else  decision  trees.  Rather  than 
representing  its  control  structure  explicitly,  AX  represents  much  of  its  control  structure  in 
terms  of  rules.  A  rule  is  a  chunk  of  knowledge  (frame)  that  has  two  parts  —  a  condition  and  an 
action.  If  the  rule's  condition  is  true,  then  its  action  can  be  executed.  For  example,  the  follow¬ 
ing  rule  suggests  when  a  P -arrival  might  really  be  a  Pg-arrival. 

RULE-1 

CONDITION: 

A  P  arrival  is  less  than  1°  from  an  event 

ACTION: 

Hypothesize  that  its  phase  should  be  Pg. 

hi  carrying  out  the  action  part  of  a  rule,  other  rules  may  be  executed.  More  than  one  rule  may 
have  its  condition  satisfied  at  any  one  time,  and  the  order  in  which  rules  are  executed  is  con¬ 
trolled  by  a  rule  interpreter.  Related  sets  of  rules  may  be  organized  into  rule-domains  that 
compartmentalize  knowledge  about  different  aspects  of  the  problem.  Each  rule-domain  can  have 
its  own  rule  format  and  its  own  style  of  rule  interpreter. 

The  advantage  of  using  rules,  rather  than  explicit  control  structures,  in  complicated  systems 
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is  that  rules  are  modular  and  relatively  independent.  Additional  capability  can  be  added  by 
adding  new  rules  to  the  appropriate  rule-domains,  with  little  regard  for  what  the  rest  of  the  sys¬ 
tem  does. 

AX  currently  has  about  20  rules  for  organizing  incoming  arrival  and  event  information  and 
hypothesizing  and  verifying  events.  About  30  other  rules  have  been  designed  for  reasoning  about 
related  arrival  groups,  hypothetical  arrivals,  and  arrival-time  residuals.  To  give  an  idea  of 
how  a  rule-based  approach  works,  the  following  section  discusses  reasoning  about  residuals. 


4.  Reasoning  About  Residuals 

During  the  association  process,  arrival  residuals  provide  evidence  for  the  existence  of 
events  and  for  deciding  which  arrivals  belong  to  which  events.  A  fair  amount  of  seismological 
knowledge  must  be  used  to  reason  about  residuals.  Although  theoretical  travel-time  curves  are 
available  for  only  certain  distance  and  depth  ranges,  arrivals  are  observed  outside  these  ranges. 
Seismic  knowledge  must  be  used  to  decide  if  a  travel-time  curve  can  be  extended.  Also,  arrivals 
have  both  arrival-time  errors  and  phase-identification  errors  that  can  often  be  identified  by 
seismologists.  As  an  example  of  how  AX  might  accomplish  similar  reasoning,  consider  the 
following  situation:  there  is  an  event  and  a  P-arrival  from  a  station  that  could  be  associated 
with  the  event,  but  it  is  beyond  the  distance  range  for  which  P  can  be  computed.  There  are 
several  possibilities  to  consider: 

(a)  The  event  has  not  been  located  well  enough. 

(b)  The  arrival  could  belong  to  another  event. 

(c)  The  arrival  could  really  be  a  PKP  or  a  diffracted  P  (PDIF)  arrival. 

(d)  The  arrival  time  could  be  wrong. 

The  following  set  of  rules  could  be  used  to  consider  whether  the  arrival  could  be  PDIF  or  not: 
RULE-2 

CONDITION: 

1)  Phase  of  arrival  is  P,  and 

2)  arrival  doesn't  match  event,  and 

3)  arrival  could  be  PDIF  from  event,  and 

4)  the  event  is  large 

ACTION: 

Hypothesize  phase  of  arrival  is  PDIF. 

RULE-3 

CONDITION: 

Magnitude  of  event  is  >  5.5 

ACTION: 

Event  is  large. 

RULE-4 

CONDITION: 

1)  Event  has  arrivals  from  more  than 
15  stations,  and 

2)  At  least  2  of  these  stations  are  at 
least  50  degrees  away 

ACTION: 

Event  is  large. 

RULE- 5 

CONDITION: 

1)  Distance  of  station  of  arrivals 
to  event  is  in  PDIF  range,  and 

2)  Travel  time  matches  PDIF 

ACTION: 

Arrival  could  be  PDIF  from  event. 

These  rules  form  the  structure  as  shown  in  Fig.  11-17. 
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The  first  rule,  RULE- 2,  states  the  conditions  necessary  to  hypothesize  that  a  P-arrival  is 
really  a  PDIF.  The  other  rules  show  how  the  subgoals  of  RULE- 2  can  be  accomplished.  RULE- 2 
uses  more  than  just  the  size  of  the  residual  to  decide  if  the  arrival  could  be  PDIF,  it  also  checks 
that  the  event  is  large  enough  to  make  observing  a  PDIF  likely.  There  are  two  ways  to  decide 
if  the  event  is  large  enough.  This  allows  the  rule  to  work  even  if  there  are  no  amplitude  data 
for  the  event. 

Once  we  have  hypothesized  that  the  arrival  is  PDIF,  we  must  still  verify  that  this  is  true, 
which  is  not  easy  since  we  have  not  ruled  out  the  other  possibilities  mentioned  above.  This  is 
handled  by  creating  a  hypothetical  arrival  frame  for  this  hypothesis: 

FRAME 

h-arrival-32 


SLOT  FACET 

AKO  value 

phase  value 

previous 


DATUM  COMMENT 

arrival- 16 

pdif  (  why  reason- 52) 

P 


The  hypothesis  derives  most  of  its  information  from  the  arrival  it  was  derived  from.  The  phase 
slot  is  tagged  with  a  comment  that  contains  the  name  of  a  frame  describing  the  reasoning  behind 
the  hypothesis. 

Once  a  hypothetical  arrival  is  created,  it  is  treated  by  the  rest  of  the  system  as  a  normal 
arrival.  Both  the  original  and  the  hypothetical  arrivals  exist  together,  each  treated  indepen¬ 
dently  by  the  system.  The  correct  decision  is  made  during  the  final  event  verification  stage. 
Several  rules  are  used  to  decide  which  arrival  is  correct.  If  the  system  can't  decide,  an  analyst 
is  asked  to  make  the  decision.  It  could  even  explain  the  reasoning  behind  the  hypothesis. 

This  type  of  reasoning  goes  well  beyond  existing  AA  programs.  They  either  make  no  adjust¬ 
ments  to  phase  identification,  in  which  case  the  arrival  may  remain  unassociated  forever,  or 
they  make  adjustments  based  only  on  residuals  from  one  event.  Nothing  is  done  to  verify  the 

adjustment.  . 

K.  R.  Anderson 
S.  T.  Rosenberg 
D.  Lanan 


E.  AN  ENTROPY  SCRAMBLED  HASH  FUNCTION 

Hash  functions  that  require  only  one  probe  per  lookup  and  a  minimum -sized  table  are  called 
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minimum  perfect  hash  functions.  ’  Although  minimum  perfect  hash  functions  can  be  designed, 
they  require  a  static  set  of  keys,  and  only  work  for  small  sets.  Thus,  their  application  is  limited. 

This  report  describes  a  new  method  of  computing  a  hash  function  from  a  subset  of  keys  that 
produces  evenly  distributed  hash  addresses  based  on  the  entropy  of  individual  bit  positions  of  the 
keys.  Although  the  method  does  not  produce  a  perfect  hash  function,  primary  clustering  is  sig¬ 
nificantly  reduced.  Thus  it  can  improve  any  hashing  method,  and  it  is  attractive  for  retrieval 
from  secondary  storage. 

i.  Entropy 

As  the  concept  of  entropy2^  is  important  in  the  analysis  of  the  hash  functions  presented  here, 
it  will  be  briefly  described.  Entropy  is  a  measure  of  information  content.  Suppose  we  had  a 
source  X  of  discrete  messages  x.  (i  =  1,...,N),  each  with  an  independent  probability  of  occurring 
p(x.).  Then,  the  entropy  of  the  information  source  X  is: 
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N 

H=—  D  p(x.)  log2[p(x.)]  bits  . 
i=l 

This  is  the  average  number  of  bits  per  message  needed  to  encode  the  messages  from  X. 

Entropy  can  also  be  thought  of  as  a  measure  of  the  shape  of  the  function  p(x).  The  more 
uniformly  distributed  p(x)  is,  the  higher  the  entropy. 

For  example,  if  X  produces  only  two  messages,  0  (with  probability  pg)  and  1  (with  prob¬ 
ability  pj  =  1  —  p0),  then  the  entropy  is: 

H  =  — pQ  log2(p0)  —  (i  —  Pq)  log-,(i  —  pQ)  .  (II-3) 

When  pQ  is  very  different  from  p^,  then  H  is  near  zero.  When  pg  is  approximately  equal  to  pj, 
then  H  is  near  its  maximum  value  of  1. 

Entropy  can  be  used  in  several  ways  to  measure  the  performance  of  hash  functions.  For 
example,  given  a  hash  function  h(k)  that  returns  an  address  in  a  table  given  a  key  k,  its  address 
entropy  is  given  by 


Ha  =-^  P(h(k)]  log2{p[h(k)]} 
k 

where  p[h(k)]  is  the  fraction  of  the  keys  that  hash  into  address  h(k).  A  perfect  hash  function 

would  distribute  its  keys  uniformly  with  one  key  per  address.  Then,  p[h(k)  ]  =  l/K,  where  K  is 

the  number  of  keys,  and 

Ha=  ^°g2(K>  • 

k 

If  a  hash  function  extracts  less  than  log2(K)  bits  of  information  from  the  keys,  then  it  will  not 

be  able  to  distinguish  between  certain  keys,  and  some  keys  will  have  the  same  initial  hash 

address. 


2.  Hash  Functions 


As  our  attention  focuses  on  the  initial  hash  address,  and  not  on  methods  of  resolving  colli- 
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sions,  a  simple  hashing  algorithm  is  considered.  It  is  similar  to  algorithm  L  of  Knuth  and 
looks  something  like  this  when  written  in  C  (see  Ref.  28); 

search(k,  keys) 

{ 

for(i  =  h(k);  keys[i]  '.=  EMPTY;  i  =  i  +  1) 
if  (keys[i)  ==  k)  return(i); 
return(  FAIL); 

} 

Given  a  key  k  and  a  table  of  keys,  keys[]  ,  search(k,  keys)  returns  the  address  in  the  table  with 
a  key  corresponding  to  k.  The  hash  function  h(k)  returns  an  address  in  the  table  given  a  key. 

The  algorithm  searches  consecutive  table  elements  starting  at  the  hash  address.  A  similar  algo¬ 
rithm  can  be  used  to  insert  elements  into  the  table.  Although  this  method  is  effected  by  primary 
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clustering,  it  is  convenient  for  searching  secondary  storage.  However,  to  make  this  method 

effective,  h(k)  must  be  designed  to  minimize  primary  clustering. 

To  see  the  effect  of  primary  clustering,  consider  a  commonly  used  hash  function  that  treats 
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the  key  as  an  integer  and  uses  it  as  the  table  address  modulo  the  table  size 
Hl(k) 

^  return(k%TABLESIZE), 

The  first  line  of  Table  II- 3  summarizes  the  results  of  retrieving  the  3553  seismic  station  codes 
from  a  table  of  size  8191  using  the  search  method  described  above.  The  address  entropy  of  Hi 
is  relatively  high  since  the  maximum  is  log2(3553)  =  11.79.  However,  HI  is  strongly  affected 
by  primary  clustering.  Although  the  table  is  only  43-percent  full,  on  the  average  it  takes  over 
10  comparisons  to  find  a  match,  and  in  the  worst  case  427  comparisons  are  made.  This  is 
significantly  worse  than  the  number  of  comparisons  predicted  theoretically  by  Knuth,  which 
is  1.37.  He  shows  that  given  the  fraction  of  the  table  that  is  full  A,  the  number  of  comparisons 
in  a  successful  search  C  is 

C  =  l+W-A)!  .  (II-4) 

The  theory  assumes  that  the  hash  function  distributes  the  keys  evenly  over  the  table.  Clearly, 
this  is  not  true  here. 


TABLE  11-3 

COMPARISON  OF  HASH  FUNCTIONS 

Table 

Unique 

Address 

No.  Comparisons 

Method 

Size 

Addresses 

Entropy 

Average 

Worst  Case 

HI 

8191 

2617 

11.23 

10.49 

427 

H2 

8192 

2346 

11.03 

2.17 

47 

H3 

8191 

2766 

11.32 

1.65 

27 

H4 

8191 

2918 

11.42 

1.35 

10 

H4 

4093 

2359 

11.04 

3.61 

112 

This  can  be  avoided  by  using  the  entropy  of  the  individual  bits  of  the  keys  to  produce  a  more 
even  distribution  of  hash  addresses.  The  bit  entropy  H.  gives  the  amount  of  information  gained 
about  a  key  by  knowing  its  i1*1  bit.  It  can  be  determined  from  Eq.  (II- 3)  where  pQ  is  the  percent¬ 
age  of  the  keys  with  that  bit  set  to  zero.  Table  II-4  shows  the  bit  entropies  for  the  bit  positions 
in  the  4-byte  station  codes.  Figure  11-18  shows  the  same  information  in  histogram  form.  From 
Table  II- 4  and  Fig.  11-18,  several  things  can  be  observed: 

(a)  Several  bits  provide  no  information  about  the  keys.  Since  the  ASCII 
character  code  is  used,  bit  7  of  each  character  is  always  0.  Also, 
since  the  first  two  characters  of  the  codes  are  always  capital  letters, 
their  sixth  bit  is  always  1. 
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(b)  Most  of  the  information  about  the  keys  comes  from  the  low-order  3  or 
4  bits  of  the  first  three  characters. 

(c)  Almost  no  information  is  contained  in  the  fourth  character.  This  is 
because  many  of  the  codes  are  only  three  letters  long. 

This  suggests  that  a  reasonable  hash  function  could  be  made  by  taking  only  the  first  1 3  bits 
of  Table  II-4, 

H2(k) 

{ 

h  =  0; 

for  (i  =  1;  i  <=  13;  i  =  i  +  1) 
h  =  h  «1  !  bit(i,  k); 
return(h); 


where  bit(i,  k)  returns  the  bit  of  the  key  k  mentioned  in  the  i4*1  line  of  Table  II-4.  The  effect  of 
the  for-loop  is  to  rearrange  the  bits  of  the  key  so  that  the  bits  with  the  highest  bit  entropy  appear 
at  the  high-order  end  of  the  key.  The  resulting  scrambled  keys  are  more  evenly  distributed 
over  the  address  space  than  are  the  original  keys. 

The  results  of  this  hash  function  are  summarized  in  the  second  line  of  Table  II- 3,  The  ad¬ 
dress  entropy  is  not  as  good  as  that  of  hash  function  HI,  because  only  13  bits  of  information  are 
used  to  compute  the  hash  function.  These  13  bits  contain  only  11.025  bits  of  independent  infor¬ 
mation  about  the  keys.  On  the  other  hand,  the  primary  clustering  has  been  significantly  reduced. 

An  obvious  improvement  can  be  made  to  the  hash  function  by  using  all  the  bits  of  a  key  that 
have  any  information  about  the  keys,  but  scrambled  by  their  bit  entropy: 

H3(k) 

{ 

h  =  0; 

for  (i  =  1;  i  <=  26;  i  =  i  +  1) 
h  =  h  «  1  !  bit(i,  k), 
return(h%TABLESIZE), 


The  corresponding  line  in  Table  II- 3  shows  that  we  get  more  than  5-percent  increase  in  the  num¬ 
ber  of  unique  addresses  over  HI  and  a  corresponding  increase  in  the  address  entropy.  The  aver¬ 
age  number  of  comparisons  and  the  worst-case  number  of  comparisons  have  been  significantly 
reduced. 

A  final  improvement  provides  an  11 -percent  improvement,  over  HI,  in  the  number  of  unique 
hash  addresses,  and  an  average  number  of  comparisons  that  agrees  well  with  Eq.  (II-4): 

H4(k) 

{ 

h  =  0; 

right  s  1.0; 

left  =  0.0; 

for  (i  =  0;  i  <=  26;  i  =  i  +  1)  { 
b  =  bit(i,  k); 

if  (b  ==  0)  right  =  left  +  (right-left)  *p0[i]; 
else  left  =  left  +  (right  -left)*(1.0  -  p0[i)); 

h  =  (left +right)*0, 5*67108864; 

return(h%TABLESIZE); 
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where  pO(i|  is  the  probability  that  bit  i  is  zero,  taken  from  Table  II-4,  and  67108864  is  2^t 
Hash  function  H4  is  similar  to  a  binary  search.  The  position  of  the  key  in  the  table  is  bounded 
by  LEFT  and  RIGHT.  If  p0[i]  was  allowed  to  be  a  function  of  the  other  bits  in  the  key,  a  binary 
search  would  result  and  a  unique  hash  address  would  be  found  for  every  key.  However,  this 
would  replace  a  hash  search  with  a  binary  search. 

The  last  line  of  Table  II- 3  shows  the  effect  of  reducing  the  table  size  in  half  on  hash  function 
H4.  The  total  number  of  unique  addresses  is  still  80  percent  of  the  value  of  the  larger  table,  and 
the  average  number  of  comparisons  agrees  well  with  that  predicted  by  Eq.  ( II - 4) .  Even  with  over 
80  percent  of  the  table  full,  its  performance  is  almost  three  times  better  than  hash  function  HI. 

In  conclusion,  the  entropy  concept  provides  several  measures  of  hash  function  performance. 
Using  bit  entropy  to  scramble  the  bits  can  provide  hash  functions  that  have  well -distributed  hash 
addresses.  The  method  is  not  limited  to  static  sets.  Usually,  the  application  can  provide  some 
information  about  how  to  scramble  the  bits.  For  example,  for  most  programming  applications, 
ordering  the  bits  so  that  the  low-order  bits  of  each  character  appear  in  the  high-order  bits  of 
the  scrambled  key  should  be  better  than  an  unscrambled  hash  function  such  as  HI. 

K.  R.  Anderson 
K.  A.  Huber 


F.  THREE-COMPONENT  POWER -DETECTION  STUDY 

Three -component  power  detection  has  been  tested  on  24  h  of  data  from  seven  stations.  Data 
used  in  the  study  were  from  2  October  1980,  which  was  the  second  day  of  the  IESD  experiment. 
The  detector  used  was  our  in-house  Z-statistic  detector  which  has  been  described  in  detail  in 
the  previous  SATS.2  Information  pertaining  to  the  stations  is  presented  in  Table  II- 5. 


TABLE  11-5 

data  AVAILABLE  FOR  THIS  STUDY 

Code 

Type 

Location 

Outages 

ANMO 

SRO 

Albuquerque,  New  Mexico 

Negligible 

BCAO 

SRO 

Bangui,  Central  African  Republic 

12.5  h  starting  at  6:00 

BOCO 

SRO 

Bogota,  Colombia 

Vanous  small  dropouts 

CP  O 

RSTN 

Cumberland  Plateau,  Tennessee 

Negligible 

HN- 

SDCS 

Houlton,  Maine 

I  h  starting  15:00 

PI- 

SDCS 

Pinedole,  Wyoming 

Various  small  dropouts 

PMR 

USSN 

Palmer,  Alaska 

Various  small  dropouts 

Analyst  picks  from  different  sources  were  available  for  the  data.  Some  of  the  analysts  re¬ 
ported  only  teleseismic  arrivals,  others  reported  everything.  Most  of  the  picking  was  done  on 
analog  records.  No  records  were  available  for  PI-,  so  the  picks  from  nearby  BDW  were  used; 
also,  no  picks  for  HN-  were  available,  so  23  h  of  digital  data  were  picked  by  hand.  Our  AA  algo¬ 
rithm  produced  23  preliminary  events  which,  in  turn,  were  used  to  generate  potential  arrival  times 
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at  the  stations.  Analyst  picks,  along  with  potential  arrival  times  generated  from  our  PEL,  pro¬ 
vided  a  means  to  test  the  detector.  Detection  logs  from  the  three  components  were  merged  and 
tested  against  the  Z -component  alone.  Table  II- 6  shows  the  results. 


TABLE  11-6 

EXPERIMENTAL  RESULTS 

Extra  H 

Analyst 

Z  Det 

3  Det 

PEL 

Z  Det 

3  Det 

Station 

Z  Dets 

Dets 

Picks 

Picks 

Picks 

Assc 

Assc 

Assc 

ANMO 

24 

9 

42 

31 

33 

23 

9 

BCAO 

10 

11 

11 

6 

6 

23 

1  I 

6 

BOCO 

5 

11 

6 

3 

5 

23 

mm 

1 

CPO 

83 

35 

2 

1 

1 

23 

m 

HN- 

5 

10 

7 

6 

6 

23 

PI- 

5) 

28 

81 

55 

58 

23 

PMR 

58 

35 

35 

34 

35 

23 

PM 

The  second  column  (Z  Dets)  in  Table  II- 6  is  the  number  of  times  the  detector  triggered  on 
the  Z -component,  the  third  column  is  the  number  of  unique  arrivals  added  by  the  horizontals. 
Clearly,  adding  H-components  can  increase  the  total  number  of  arrivals  anywhere  from  37  to 
200  percent.  Column  4  contains  the  number  of  analyst  picks,  while  column  5  contains  the  num¬ 
ber  of  those  picks  that  were  within  the  detector  on  and  off  times  on  the  Z -component.  The  sixth 
column  contains  similar  information  for  all  three  components.  At  three  of  the  stations,  adding 
horizontals  made  no  difference  as  far  as  the  number  of  analyst  picks  detected,  at  the  other  four, 
the  difference  varied  from  one  more  detected  at  PMR  to  three  more  detected  at  PI-.  Column  7 
shows  the  number  of  potential  arrivals  generated  from  the  PEL.  The  number  of  associated  ar¬ 
rivals  detected  on  the  Z-component  are  shown  in  the  eighth  column,  while  the  ninth  column  shows 
the  corresponding  information  for  all  three  components.  Adding  the  detections  from  the  hori¬ 
zontals  increased  the  number  of  associated  arrivals  detected  at  four  of  the  stations. 

Certainly  3-component  detection  seems  to  provide  us  with  more  significant  arrivals,  but 
unfortunately  there  is  also  a  three-fold  increase  in  false  alarms.  One  would  expect  this  from 
a  Gaussian  distributed  detection  statistic  (log  of  the  average  power).  Figure  11-19  shows  the 
distribution  of  the  Z-statistic  (0  mean,  S. D.  =  1). 

Figure  II-20(a-c)  shows  some  unassociated  arrivals  that  were  only  detected  on  the 
H-components.  Figure  II-20(a)  shows  a  short- period  Love  wave. 

Figures  11-21  to  11-23  show  the  associated  waveforms  for  three  events  which  had  some  asso¬ 
ciations  made  from  H-component  detections.  Figure  II-21(a-b)  shows  event  1108,  which  is  a 
split-off  from  a  Fiji  event.  The  split  would  have  occurred  with  or  without  the  H-component 
detections.  Both  the  ANMO  and  PMR  arrivals  were  detected  on  the  horizontals  alone.  A  split- 
off  from  an  Eastern  USSR  event  also  occurred;  however,  this  split  could  not  have  been  formed 
without  the  horizontals  (this  event  is  not  shown).  Figure  II-22(a-b)  shows  event  1109,  located 


in  Colorado,  which  was  not  reported  by  USGS  or  Sweden.  Without  the  detections  on  the  horizon¬ 
tals  at  ANMO,  there  would  have  been  insufficient  arrivals  to  associate  or  locate  this  event. 
Event  1113,  located  in  Western  Idaho,  is  shown  in  Fig.  II-23(a-b);  this  event  is  another  original 
produced  during  the  first  pass  of  our  AA.  If  not  for  the  detections  on  the  horizontals  at  BCAO, 
this  event  could  not  have  been  associated  or  located.  It  is  interesting  to  note  that  the  detector 
missed  the  detection  on  the  Z- component,  while  it  did  manage  to  detect  it  on  the  horizontals. 

Thus,  we  find  that  adding  the  detections  from  the  horizontals  can  increase  the  number  of 
significant  arrivals  for  a  given  period  of  time.  And,  while  one  can  argue  that  adding  two  mar¬ 
ginal  events  and  one  split  event  to  an  event  list  is  questionable  given  the  requirement  to  process 
three  times  as  many  false  alarms,  it  seems  that  with  additional  refinement  3-component  detec¬ 
tion  can  become  a  significant  enhancement  to  the  Automatic  Signal  Processor. 

M.  A.  Tiberio 


C..  COMPARISON  OF  ANALYST  AND  AUTOMATIC  TECHNIQUES 

FOR  PICKING  FIRST  ARRIVALS 

As  part  of  the  U.S.  contribution  to  the  International  Data  Collection  Experiment,  several 
methods  of  picking  first  arrivals  were  compared  for  the  number  of  associated  arrivals  generated 
and  accuracy  of  picks.  The  database  used  consists  of  all  the  digital  and  analog  data  collected 
on  2  and  4  October  1980  from  the  Global  Digital  Seismograph  Network  (GDSN),  as  well  as  re¬ 
ported  readings  from  selected  WWSSN  stations.  The  NEIS  monthly  summary  provides  a  check 
of  each  method  by  identifying  which  generated  arrivals  are  associated  with  known  events.  The 
criteria  for  association  of  an  arrival  with  a  given  event  are  very  broad  —  an  arrival  is  considered 
associated  if  its  time  agrees  to  within  1  min.  with  that  predicted  by  the  Jeffreys -Bullen  travel¬ 
time  tables.  Such  a  range  of  acceptability  is,  of  course,  far  from  being  tolerable  for  event  loca¬ 
tion,  but  was  found  necessary  to  associate  arrivals  from  particular  events,  which  for  some  sta¬ 
tions  had  large  residuals  no  matter  which  method  was  used  to  make  the  arrival-time  pick.  The 
difference  between  the  time  pick  and  the  arrival  time  predicted  by  the  Jeffreys -Bullen  time  tables 
provides  a  measure  of  the  accuracy  of  each  method. 

The  four  methods  compared,  ranging  from  completely  manual  to  fully  automatic  detection 
and  analysis,  are  listed  below. 

(1)  Analyst  ASL  Detector  —  Analyst  review  of  digital  detections  made  by  the 
Albuquerque  Seismological  Laboratory  (ASL)  or  SRO  detector.  The  dig¬ 
ital  data  are  plotted  on  a  large  scale,  and  time  and  amplitude  measure¬ 
ments  made  by  very  precisely  using  a  hardware  cursor. 

(2)  Automatic  LL  Detector  —  Results  of  automatic  processing  as  described  in 
Sec.  I-D  of  this  SATS.  In  cases  where  the  automatic  time  picker  failed 
even  though  the  "Z"  -detector  had  identified  a  signal  which  was  in  fact 
associated,  the  time  is  that  at  which  the  detector  triggered.  No  analyst 
review  or  interaction. 

(3)  Analyst  SRO  Analog  -  Analyst  scanning  film  records  of  SRO/ASRO  data, 
performing  detection  and  measurements  manually. 

(4)  Analyst  Station  (ALQ)  —  Analyst  measurements  made  as  (3)  but  on  an  alter¬ 
nate  station  (e.g.,  ALQ  WWSSN  station  as  alternative  to  ANMO).  A  con¬ 
ceivable  advantage  to  such  picks  is  that  they  are  made  at  the  station  itself 
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by  an  analyst  with  considerable  experience  in  recognizing  signals 
particular  to  both  the  site  and  specific  source  regions. 

One  advantage  of  the  automatic  processor  is  in  the  speed  with  which  the  data  can  be  analyzed. 
Detection,  windowing,  timing,  and  amplitude/period  measurement  on  24  h  of  continuous  short- 
period  (20  samples/s)  data  requires  less  than  15  min.  of  computer  time  on  a  moderate-sized 
computer.  This  is  considerably  faster  than  a  human  analyst  processing  the  same  number  of 
arrivals  (~30/day)  from  analog  records. 

Techniques  (1)  through  (4)  above  generated  98,  111,  89,  and  69  associated  arrivals,  respec¬ 
tively.  Of  these,  66,  70,  63,  and  52  (67,  63,  71,  and  75  percent),  respectively,  had  travel-time 
residuals  of  less  than  2  s.  The  distributions  of  the  residuals  by  each  technique  are  given  in 
Table  (1-7.  It  can  be  seen  that  none  of  the  four  methods  appears  to  produce  picks  which  are 
either  consistently  late  or  early.  The  automatic  technique  is  slightly  more  prone  to  generating 
picks  with  residuals  greater  than  10  s,  which  is,  of  course,  unacceptable  for  the  purposes  of 
accurate  location. 

Moderate  travel-time  residuals  (e.g.,  <3  to  4  s)  do  not  necessarily  indicate  a  poor  arrival¬ 
time  pick  but  can  be  due  to  regional  variations  in  the  travel-time  pick.  In  order  to  attempt  to 
cancel  out  such  path-specific  bias,  the  difference  between  the  times  determined  by  methods 
(2)  through  (4)  above  and  that  produced  by  (1)  has  been  calculated  for  the  associated  picks  only. 


TABLE  11-7 

DISTRIBUTION  OF  TRAVEL-TIME  RESIDUALS 

FOR  ASSOCIATED  ARRIVALS  BY  EACH  METHOD 

Residual  R 

Analyst 

Automatic 

Analyst 

Analyst 

(s) 

ASL  Detector 

LL  Detector 

SRO  Analog 

ST  A  Analog 

R  <  -20 

i 

-20  <  R  <  -10 

2 

2 

i 

-10  <  R  <  -5 

1 

' 

2 

1 

-5  <  R  <  -4 

3 

0 

3 

0 

-4  <  R  <  -3 

2 

2 

2 

0 

-3  <  R  <  -2 

6 

9 

3 

3 

-2  <  R  <  -1 

12 

16 

9 

8 

-1  <  R  <  0 

22 

29 

19 

17 

0  <  R  <  1 

20 

13 

28 

15 

1  <  R  <  2 

12 

12 

7 

12 

2  <  R  <  3 

8 

8 

5 

2 

3  <  R  <  4 

2 

3 

1 

1 

4  <  R  <  5 

2 

4 

3 

2 

5  <  R  <  10 

3 

3 

4 

5 

10  <  R  <  20 

1 

5 

2 

2 

R  >  20 

2 

3 

0 

1 
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TABLE  11-8 

DISTRIBUTION  OF  DIFFERENCES  BETWEEN  TIME  PICKS  MADE 

BY  METHODS  (2),  (3),  AND  (4)  AND  METHOD  (1) 

Difference  D 

Methods 

« 

- 

(2)  and  (1 ) 

(3)  and  (1 ) 

(4)  and  (1) 

D  <  -10 

2 

1 

0 

-10  <  D  <  -5 

2 

0 

0 

-5  <  D  <  -2 

3 

1 

0 

-2  <  D  <  -1.5 

2 

1 

1 

-1.5  <  D  <  -1.0 

6 

2 

1 

-  1.0  <  D  <  -0.5 

12 

2 

3 

-0,5  <  D  <  0 

22 

12 

4 

0  <  D  <  0.5 

26 

15 

21 

0.5  <  D  <  1.0 

1 

2 

4 

1.0  <  D  <  1.5 

3 

1 

1 

1.5  <  D  <  2.0 

3 

0 

0 

2.5  <  D  <  5.0 

4 

0 

3 

5.0  <  D  <  10.0 

1 

1 

1 

D  >  10.0 

4 

1 

2 

(1)  =  Analyst  pick,  ASL  detection 

(2)  =  Automatic  pick,  LL  detection 

(3)  =  Analyst  pick,  SRO 

analog 

(4)  =  Analyst  pick,  alternate  station  (9  stations  only  of  total  17) 

Technique  (1)  was  chosen  as  the  standard  since  it  is  reasonable  to  assume  that  the  ability  to 
plot  the  seismogram  on  a  very  large  scale  from  the  digital  data  will  enable  the  most-accurate 
onset  time  picks  to  be  made  by  the  analyst.  The  results  are  shown  in  Table  II-8,  The  tendency 
of  automatic  technique  (2)  to  produce  more  erroneous  picks  than  (3)  or  (4)  is  still  apparent,  as 
is  its  ability  to  produce  more  total  picks  than  method  (3)  which  uses  analog  data  from  the  same 
seismometer.  Perhaps  the  most  surprising  result  is  the  occasional  large  difference  between 
the  times  picked  by  the  various  human  analysts.  Even  for  arrivals  designated  by  various  ana¬ 
lysts  as  impulsive,  denoting  an  accuracy  which  the  analyst  believed  to  be  +0.2  s,  differences  of 
1  s  or  larger  are  quite  common. 

The  GDSN  studied  here  consisted  of  17  stations  at  the  time  of  the  data-collection  experiment. 
For  the  32  events  on  the  monthly  summary,  9  were  recorded  to  minimum  distances  of  less  than 
10*.  Of  the  remaining  23,  only  10  were  detected  at  5  or  more  stations  of  the  GDSN  network  and 
could  conceivably  be  located  by  data  from  these  stations  if  the  azimuthal  coverage  were  adequate. 
Addition  of  data  from  the  7  other  stations  within  the  Continental  U.S.  and  Alaska  may  improve 
this  situation  slightly.  R.  G>  North 

M.  W.  Shields 
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Fig.  II- 3.  Computer  description:  50.55  s  of  noise,  a  glitch  0.1  s  wide 
and  —  560  units  high,  10.1  s  of  noise,  a  box  with  a  height  of —  1091  be¬ 
ginning  with  a  glitch,  and  115.25  s  of  noise. 
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Fig.  II-4.  Detail  of  a  calibration  pulse  (solid  line)  and  its  hysteresis 
filtered  version  (dotted  line). 
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Fig.  II- 7.  Dynamic  match  between  two  dissimilar  waveforms. 
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Fig.  II- 8.  Cluster  analysis  of  envelopes  of  eight  local  seismograms 
recorded  at  BDW.  A  dendrogram  indicating  hierarchical  clustering 
of  waveforms  is  shown  by  dotted  lines  to  left  of  envelopes. 
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Fig.  II- 9.  Tree  structure  of  affinity  based  on  average  magnitude  criteria. 


Fig.  II' 10.  Representation  of  waveform 
at  different  levels  down  tree  based  on 
average  amplitude  criteria:  (a)  level  3, 
(b)  level  5,  and  (c)  level  7. 
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Fig.  11-15.  Inheritance  hierarchy 
of  frames  representing  seismic 
parameters. 


Fig.  11-16.  Inheritance  hierarchy 
of  frames  representing  seismic 
waves. 
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Fig.  11-17.  Structure  formed  by  set  of  five  rules  used  to  determine 
whether  or  not  an  arrival  could  be  PDIF. 
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Fig.  II-  18.  Bar  graph  showing 
relative  entropies  of  each  bit 
of  4-byte  keys.  Height  of  bars 
corresponds  to  divisions  of 
Table  II-4. 
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STANDARD  DEVIATIONS 


Fig.  II- 19.  Z -statistic  distribution. 


pmr 
izrsu 
17010  960 
2000  nm 


pmr 
smsu 
17010  960 
2700  nm 


pmr 
sersu 
17010.960 
49.00  nm 


hn- 

izrsn 

60494  950 
33.00  nm 


hn- 

$nrsn 

60494.950 
119  00  nm 


hn- 

*«r*n 

60494  950 
153  00  nm 


5! 


(C) 

Fig.  II-20(a-c).  Continued. 
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Fig.  Il-21(a-b).  Bogus  Alaskan  event  split  from  Fiji  event. 
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Pig.  II-22(a-b).  Hypothetical  Colorado  event. 
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III.  GENERAL  SEISMOLOGY 


A.  A  REVIEW  OF  FREQUENCY- DEPENDENT  ANELASTICIT Y 

1.  Frequency- Dependent  Anelasticity  from  Spectral  Measurements 

Many  studies  have  shown  that  the  frequency  content  of  body  waves  in  the  band  0.5  to  10  Hz 
is  too  high  to  be  consistent  with  the  intrinsic  anelasticities  and  determined  from  observa¬ 
tions  of  seismic  waves  in  the  band  0.005  to  0.1  Hz  (see  Refs.  1  and  2).  The  results  of  these  stud¬ 
ies  are  confirmed  by  the  spectra  of  short-period  P-waves  observed  at  SRO/ASRO  stations  from 
deep  focus  events.  The  body-wave  spectra  of  deep  events  are  desirable  in  attenuation  studies 
because  they  are  uncontaminated  by  near- surface  multiples,  and  they  do  not  include  the  effects 
of  anomalous  elastic  and  anelastic  structures  that  may  be  present  in  the  source  region.  Eight 
events  were  chosen  for  study  having  foci  from  350  to  650  km  depth  and  body-wave  magnitudes 
from  5.5  to  m^  6.0.  Spectral  studies  by  Wyss  and  Molnar3  of  deep  focus  events  of  this  size 
found  the  P-wave  corner  frequency  to  range  from  0.1  to  0.6  Hz.  This  corner  frequency  is  gen¬ 
erally  less  than  the  frequency  band  of  short-period  SRO/ASRO  stations  in  which  the  signal-to- 
noise  ratio  (SNR)  is  high  for  teleseismic  body  waves.  A  common  practice  is  to  assume  an  u.’”2 
rolloff  above  the  corner  frequency  and  to  attribute  the  additional  decay  of  the  observed  spectra 
to  the  effects  of  intrinsic  anelasticity.  The  spectra  of  eight  events  were  computed  by  correcting 
for  the  instrument  response  and  applying  a  series  of  narrowband  Gaussian  filters  equivalent  to 

4 

the  MARS  analysis  described  by  Savino  et  al.  Filter  outputs  that  were  less  than  half  the  long¬ 
term  average  of  the  time-domain  output  of  the  filter  or  that  lay  outside  of  an  acceptance  window 
for  group  arrival  time  were  rejected  as  noise.  Figure  111- 1  shows  the  spectra  and  anelastic 
interpretations  computed  for  a  deep  f  jcus  event  in  Peru  recorded  at  ANMO.  For  a  spectral 
rolloff  of  o.-2,  a  t*  =  0.2  is  inferred,  (t*  is  the  path  integrated  effect  of  anelasticity  on  a  body 
wave  and  is  commonly  defined  by  t*  =  J Q-1dt,  where  the  travel-time  integral  is  taken  along  the 
ray  path.) 

Assuming  double  t*  for  an  equivalent  surface  focus  event  gives  t*  =  0.4.  Figure  III- 2  shows 

a  -2 

the  spectra  for  a  shield  station  NWAO.  Here,  the  interpretation  with  an  u>  rolloff  gives 

t*  =  0.05  (O.i  for  an  equivalent  surface  focus  event).  All  the  events  required  a  t*  ^  0.2,  con- 
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firming  the  results  of  Der  et  al.  that  anelasticity  must  increase  with  frequency  in  order  to 

account  for  the  observed  spectra  of  short-period  body  waves. 

2.  ScS/ScP  Synthetics:  Implications  for  a  High-Frequency  Cutoff 
in  the  Seismic  Relaxation  Spectrum 

Burdick^  determined  the  time-domain  operator  needed  to  convert  an  ScP  phase  into  an  ScS 
phase  recorded  at  the  same  receiver  site,  and  inferred  from  these  results  the  high-frequency 
limit  i/2irTm  of  the  relaxation  spectrum  of  the  mantle.  A  fundamental  observation  to  match  in 
his  data  set  of  short-period  WWSSN  records  is  that  the  apparent  period  of  ScS  is  nearly  always 
at  least  twice  that  of  ScP.  Since  both  phases  leave  the  source  as  S-waves  at  a  takeoff  angle  that 
differs  by  less  than  9”,  the  difference  in  frequency  content  can  be  primarily  attributed  to  attenu¬ 
ation  in  the  mantle  rather  than  to  source  directivity.  Burdick  tested  a  model  in  which  the  fre¬ 
quency  dependence  was  assumed  to  be  constant  with  depth.  With  this  model  he  concluded  that 
strong  frequency  dependence  in  the  mantle  does  not  occur  until  frequency  exceeds  1  Hz.  This 
result  would  imply  that  anelastic  models  giving  t*  =  1  are  still  appropriate  for  body  waves 


observed  on  short-period  WWSSN  records  having  dominant  energy  at  1  Hz.  Using  this  result 
and  assuming  an  absorption-band  model  that  allows  Q  to  increase  linearly  with  frequency  above 
1  Hz,  it  is  impossible  to  match  the  high-frequency  content  of  the  ScS  phases  observed  by  Sipkin 
and  Jordan  and  the  body  waves  that  bottom  in  the  mid-mantle  observed  by  Der  et  al.5 

Experiments  with  synthetic  ScS  and  ScP  phases  in  both  the  time  and  frequency  domains  were 
conducted  in  order  to  confirm  Burdick's  results  and  to  investigate  the  effects  on  ScS  and  ScP  of 
source  pulse  width,  instrument  band  width,  and  the  effects  of  anelastic  models  that  vary  with 
depth  as  well  as  frequency.  Spectra  and  seismograms  were  synthesized  by  the  methods  de- 
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scribed  in  Lundquist  and  Cormier  for  models  of  the  relaxation  spectrum  of  the  mantle.  All 

calculations  were  performed  in  the  isotropic  PREM  model9  referenced  at  1  Hz,  the  attenuation 

model  of  which  gives  t*  =  1  for  a  teleseismic  P-wave  of  an  event  having  a  surface  focus.  The 

lower  time  constant  of  a  mantle-wide  absorption  band  was  varied  from  t  =  0.01  to  0.9  s/rad. 

m 

The  frequency- independent  Q  model  of  PREM  was  reproduced  at  frequencies  below  16  Hz  for 
Tm  =  and  below  0.06  Hz  for  =  0.9  s/rad.  Time-domain  pulses  were  synthesized  for  a 

600-km-deep  event  using  a  triangle  source  pulse  with  a  width  that  varied  between  0.25  to  2  s. 
Figure  ITT—  3  shows  the  results  for  ScS  and  ScP  phases  observed  with  a  WWSSN  short-period 
instrument  for  a  1-s  source  pulse.  For  all  source  pulse  widths,  it  was  impossible  to  obtain 
the  observed  shift  in  frequency  of  ScS  relative  to  ScP  without  having  Tm-s:  0.  1,  which  corre¬ 
sponds  to  a  t*  0.  7  at  1  Hz  for  a  surface  focus  P-wave  observed  at  60°.  This  test  essentially 
confirms  the  results  of  Burdick  for  a  mantle-wide  absorption  band. 

Another  series  of  tests  was  performed  for  a  mantle-wide  relaxation  band  that  possessed 
a  high  attenuation  zone  at  the  base  of  the  mantle  of  the  type  suggested  by  Mitchell  and 
Helmberger,10  which  has  =  100  and  =  300  in  a  150-km  zone  above  the  core-mantle  bound¬ 
ary.  Figure  HI-4  compares  the  results  for  the  synthetic  spectral  ratio  ScS/ScP  for  this  model 
(dotted  lines)  with  the  results  of  a  model  having  small  attenuation  at  the  base  of  the  mantle  (solid 
lines).  The  effect  of  a  low  Q  zone  of  this  magnitude  is  not  sufficient  to  move  the  models  having 
small  attenuation  at  high  frequencies  (r  >0,1)  into  the  domain  of  ScS/ScP  observations. 

A  third  type  of  relaxation  model  was  tested  in  which  the  relaxation  band  was  assumed  to  lie 
in  the  high-frequency  band  in  the  upper  400  km  and  lowermost  150  km  of  the  mantle,  but  shifts 
to  frequencies  less  than  0.06  Hz  in  the  mid-mantle.  In  the  upper  mantle,  this  model  resembles 
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the  double  absorption  band  proposed  by  Lundquist  and  Cormier  to  satisfy  the  waveforms  of 
short-period  body  waves.  The  additional  absorption  band  at  high  frequencies  at  the  base  of  the 
mantle  follows  the  model  used  by  Anderson  and  Given  to  simultaneously  satisfy  free-oscillation 
and  body-wave  data.  They  justify  the  depth  dependence  of  a  seismic  absorption  band,  in  which 
attenuation  is  shifted  to  low  frequencies  in  the  mid-mantle,  by  the  effects  that  the  pressure  and 
homologous  temperature  profiles  are  expected  to  have  on  physical  mechanisms  of  anelasticity. 
With  this  type  of  triple  absorption-band  model,  it  is  possible  to  satisfy  the  relative  frequency 
content  of  ScS  and  ScP  phases  (dashed  lines  in  Fig.  Ill- 4)  while  obtaining  a  t*  =  0.5  Hz  for  a 
surface  focus  P-wave  bottoming  in  the  mid-mantle.  A  more  intense  or  thicker  zone  of  attenua- 
tion  at  the  base  of  the  mantle  than  that  proposed  by  Mitchell  and  Helmberger,  however,  is 
required  in  order  to  bring  a  curve  having  Tm  =  0.4  into  the  domain  of  ScS/ScP  observations. 

An  experiment  similar  to  ScS/ScP  in  eliminating  the  source  spectrum  is  that  of  sS/sP. 
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For  sS,  sP  phases  bottoming  in  the  upper  700  km  of  the  mantle,  Burdick  needed  to  satisfy  an 
observed  shift  to  lower  frequency  in  the  unconverted  sS  phase.  Even  though  these  observations 
were  made  on  narrowband  instruments,  it  would  be  difficult  to  satisfy  the  observed  frequency 
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shift  without  placing  the  absorption  band  in  the  instrument  band.  These  results  demonstrate 
the  difficulty  in  simultaneously  satisfying  the  spectral  content  of  body  waves  bottoming  in  the 
mid-mantle,  and  the  ratios  ScSn+*/ScSn,  SeS/ScP,  and  sS/sP,  without  allowing  the  limits  of 
the  seismic  relaxation  band  to  strongly  vary  as  a  function  of  depth.  The  simplest  model  that 
can  satisfy  all  data  is  one  having  absorption  in  the  high-frequency  band  of  seismic  body  waves 
in  the  upper  and  lowermost  mantles,  but  shifted  to  frequencies  less  than  0.1  Hz  in  the  mid¬ 
mantle.  Regional  variation  of  ScSn+*/ScSn  observed  by  Sipkin  and  Jordan*5  and  amplitude 
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variations  in  core  grazing  PKP  waves  observed  by  Sacks  et  al.  suggest  regional  variation  in 
the  absorption  bands  in  both  the  upper  and  lowermost  mantles. 

Experiments  that  eliminate  the  source  function  can  provide  additional  constraints  on  the 
interpretation  of  body-wave  spectra  in  terms  of  absorption-band  models.  These  experiments 
can  now  be  repeated  with  greater  precision  using  digital  data  recorded  over  a  broader  band. 

The  combined  interpretation  of  source-canceling  experiments  and  spectra  of  short-period  body 
waves  should  include  the  type  of  depth  dependence  of  an  absorption  band  that  may  be  expected 
for  the  temperature  and  pressure  profile  of  the  mantle.  Quite  different  depth  dependences  of 
the  absorption  bands  can  satisfy  the  spectral  ratios  of  a  source-canceling  experiment.  It  may 
be  possible,  however,  to  show  that  for  a  given  seismic  moment  and  source  time  function  a  par¬ 
ticular  absorption-band  model  that  satisfies  spectral- ratio  data  would  produce  a  body  pulse  that 
is  too  small  to  be  observed.*5  It  is  thus  essential  for  time-domain  modeling  to  model  absolute 

as  well  as  relative  amplitudes.  ,,  „  _ 

r  V.F.  Cormier 


B.  THE  CALCULATION  OF  dA/dp  AND  OF  PARTIAL  DERIVATIVES 
FOR  TRAVEL-TIME  INVERSION  IN  TRANSVERSELY  ISOTROPIC 
SPHERICAL  EARTH  MODELS 

The  increasing  evidence  for  anisotropy  in  the  Earth's  upper  mantle  and  the  fact  that  the  re¬ 
cent  Preliminary  Reference  Earth  Model  (PREM)*^  incorporates  a  transversely  isotropic  region 
make  it  necessary  to  take  into  account  anisotropy  in  ray  theoretic  calculations,  such  as  the 
calculation  of  travel-time  curves,  geometrical  spreading  factors,  and  in  the  inversion  of  travel 
times.  The  most  general  kind  of  anisotropy  which  can  exist  in  a  spherically  symmetric  Earth 
model,  corresponding  to  the  spherically  averaged  Earth,  is  transverse  isotropy,  in  which  the 
tensor  of  elastic  parameters  is  invariant  under  rotations  about  any  radius  vector  and  is  char¬ 
acterized  by  the  five  elastic  parameters  N,  L,  A,  C,  and  F,  each  functions  of  radius  r.  It  is 
the  purpose  of  this  contribution  to  give  formulas  for  the  calculation  of  dA/dp  for  such  a  model, 
and  also  to  present  results  for  the  partial  derivatives  of  travel  times  with  respect  to  model 
parameters  which  are  needed  in  linearized  travel-time  inversion. 

Results  for  travel  time  T  and  distance  A  as  functions  of  ray  parameter  p  have  been  given 
by  Woodhouse*7  and  are  repeated  here.  The  one-way  intercept  time  may  be  written: 


T<P)  =  J  qT(P,  r)  dr 

and 

Tip)  =  J  qT(p,  r)  dr 


(III-l) 


Aip) 


(III-2) 
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where 


=  -a<lT/aP 


qT  =  qT  +  pqA 


(lit—  3 ) 


In  the  above,  and  in  later  integrals,  it  is  understood  that  the  integration  is  over  intervals  of 

radius  for  which  q  is  real,  and  that  for  rays  with  multiple  legs  the  appropriate  multiple  con- 

T  17 

tributions  to  the  integrals  should  be  taken.  The  integrands  are  given  by  : 


For  SH 


/  M  2,1/2 


Qrr  " 


-1  . 


'jb-i 


r  q 


and  for  P-SV  (upper  and  lower  signs,  respectively) 
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and 


Jr  =  (S1  “  T  R 
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'h  -  /;["! '  k  (*A  */)] 
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s  p4 
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+  2s-  X—  + 
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si  =  T(l  + 

s2  =  T(£  - 


s3  =  ~  (AC  -  F2  -  2LF) 


2  A 
s4  =  s3  -  r 


s5  -  t£(1  +  t’  "  S1S3 


(III-4) 


(in-5) 


(311-6) 


(111-7) 


where  p  is  density  (it  may  be  noted  that  the  integrands  depend  only  upon  the  ratios  L/p,  N/p, 
A/p,  C/p,  and  F/p  which  have  the  attributes  of  squared  wave  velocities). 
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1.  The  Calculation  of  dA/dp 

A  turning  radius  b(p)  for  a  ray  with  ray  parameter  p  is  determined  by 


q  (p,  r)  =0 

I r=b(p) 

and  it  will  be  assumed  that 


(II1-8) 


q  >0  for  r  >  b 


For  the  sake  of  definiteness  and  simplicity,  we  shall  restrict  our  attention  here  to  a  single 
smoothly  varying  region  of  a  spherical  model  (c  <  r  <  a)  in  which  a  turning  radius  b(p)  exists 
(c  <  b  <  a).  The  integrands  q^,  q^.  are  singular  at  r  =  b,  and  in  order  to  evaluate  the  integrals 
numerically  the  integrands  must  be  made  regular  there.  A  simple  way  to  regularize  them  is 
to  define  a  new  independent  variable 

u=(r-b)1^2  ,  i.e.,  r  =  b(p)  +  u2  .  (III-9) 

Then, 


A(p)  =  f 


(a-b) 
o 


1/2 


2G(p,  b  +  u  )  du 


(III-10) 


where 

G(p,r)  =  (r-b)l/2  qA(p,r)  (III-ll) 

and  similarly  for  Tip).  The  singularity  in  the  integrand  has  been  eliminated  and  the  integral 
may  be  evaluated  numerically  by  any  suitable  method  (e.g.,  Romberg  or  Gaussian  integration). 
To  evaluate  dA/dp,  we  write 


Aip) 


■i 


b(p)  [r  -  b(p)] 


(III— 12) 


Differentiation  under  the  integral  sign  is  not  admissible  directly  because  of  the  singularity 
at  the  turning  radius;  however,  we  may  write 


Aip)  =  f  -;(P‘-r) — dr  +  2(a-b)1/2  G(p, 
(r  — 


b) 


(III-13) 


where  the  second  term  arises  from  explicit  integration  of  the  term  subtracted  from  the  integrand. 
The  new  integrand  vanishes  at  the  turning  point,  and  differentiating  under  the  integral  and  sim¬ 
plifying  we  find: 


3A 

Tp  ‘ 


- ^72  (G  <p,  r)  +  4b'(P> 

(r  —  b)1' 2  p  2 


G(p,  r)  —  G(p,  b)  ^  b'(p)  G(P,  b) 

r“b  ”(a-b),/z 


(III-14) 


where  subscript  p  indicates  partial  differentiation.  The  integrand  here  has  an  integrable  sin¬ 
gularity  which  may  be  regularized  by  the  transformation  (III-9).  It  may  be  shown  that  the  quan¬ 
tities  b'(p)  and  G(p,b)  in  Eq.(IlI-14)  are  given  by 
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V 


b'(p)  = 


-dqyap 


8qT2/8r 


r=b 


G(p.b)  =  |-b'(p)(»q^/Sr)l/2 


r=b 


(III-15) 


It  is  straightforward  to  obtain  expressions  for  the  various  quantities  appearing  in  Eqs.  (HI-14) 
and  (III-15)  from  Eqs.  (III-ll),  (III-4),  (III- 5),  (III—  6),  and  (III— 7 );  since  they  are  lengthy,  we  will 
not  list  them  explicitly.  It  may  be  noted  that  the  derivatives  of  the  model  parameters  with  re¬ 
spect  to  radius  are  required  at  the  turning  point,  and  should  be  evaluated  (preferably  analyti¬ 
cally)  in  terms  of  the  interpolants  used  to  define  the  model.  With  a  model,  such  as  PREM, 
specified  in  terms  of  low-order  polynomials  in  radius  this  is  a  simple  matter,  but  for  a  model 
specified  pointwise  a  suitable  interpolation  scheme  (e.g.,  cubic  splines)  must  be  chosen. 


2.  Partial  Derivatives  with  Respect  to  Model  Parameters 

Let  us  suppose  that  the  Earth  model  is  specified  by  a  number  of  parameters  tj..  These  may 
be,  for  example,  polynomial  coefficients  or  values  of  wave  velocities  at  specific  radii.  Having 
specified  an  interpolation  scheme,  we  have 

A(r)  =  A(r,»].) 


C(r)  =  C(r,  rjj)  etc. 

and 

qT  =  q^p.r.rjj) 
qT  =  qT(p.  r,T|.) 

qA  =  qA(p#r*’,i)  •  (III- 16) 

In  order  to  incorporate  travel-time  data  into  a  linearized  inversion  technique  for  the  estima¬ 
tion  of  the  model  parameters  rjj,  it  is  necessary  to  find  expressions  for  (3T/8r))^,  where  the 

subscript  indicates  the  variable  to  be  held  constant  in  the  partial  differentiation  and  where  77 

1 8 

denotes  any  one  of  the  parameters  tj j.  Julian  and  Anderson  have  shown  in  the  case  of  isotropy 
that 


( — )  =  ( — ) 
18t)  A  '  8  tj  p 


(III- 17) 


where  r  is  the  intercept  time  [Eq.  (Ill— 1)],  and  this  result  can  be  immediately  extended  to  the 
present  case;  we  have 


dT  =  (|I)  dp  +  (§2J)  dr, 
dp  V  3  j)  p 


1  p  tj  8A  tj 
using  p  =  dT/dA.  Thus, 

/ 8T ,  ,8T.  .  8A,  ,8p. 

(¥^a  -  (^P  +  p(‘5p,ti  a 


(in-18) 
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But, 


so  that 


whence 
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1  ep  'r? 
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(III- 19) 


as  stated  in  Eq.  (Ill- 17). 

Using  Eq.  (ITT- 47 )  it  is  straightforward,  but  again  algebraically  involved,  to  calculate  the 
partial  derivatives  of  travel  times  with  respect  to  a  model  parameter  tj.  We  have  simply: 


(9T,  B 

{9n  a 


dr  ; 


the  integrand  is  singular  but  integrable,  and  the  integrals  may  be  evaluated  using  the  transforma¬ 
tion  Eq.  (III-9).  This  causes  no  difficulty  if  it  is  desired  simply  to  estimate  the  model  parame¬ 
ters  tj j  as  implicitly  assumed  here.  Since  the  integrands  are  not  square  integrable,  however, 

it  is  not  possible  to  apply  directly  the  apparatus  of  Backus-Gilbert  inverse  theory.  The  solution 
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to  this  difficulty  has  been  given  by  Johnson  and  Gilbert  who  have  shown  that  integration  by 
parts  may  be  used  to  recast  the  problem  in  terms  of  square  integrable  kernels  -  a  technique 
which  could  also  be  applied  in  the  present  case. 

A  computer  program  has  been  developed  to  perform  the  calculations  outlined  here,  for  an 
arbitrary  seismic  phase  and  source  depth,  under  the  assumption  that  the  model  is  parameterized 
in  the  same  way  as  the  PREM  model  of  Dziewonski  and  Anderson.1^ 

J.  H.  Woodhouse 
T.  P.  Girnius 


C.  THE  EXCITATION  OF  SEISMIC  WAVES  BY  A  SOURCE  POSSESSING 
NON- NEGLIGIBLE  SECOND  GLUT  MOMENTS 

The  technique  of  Dziewonski,  Chou,  and  Woodhouse2^’21  for  the  routine  estimation  of  seis¬ 
mic  moment  tensors,  centroid  locations,  and  centroid  times,  using  long-period  SRO  data  has 
been  extended  in  several  ways  in  order  to  gain  more  information  about  seismic  sources  not 

adequately  modeled  by  point  sources  in  space  and  time.  These  extensions,  described  here  and 
20  22 

in  previous  SATSs,  ’  have  included  the  estimation  of  source  duration  and  rupture  velocity, 
the  simultaneous  inversion  for  a  number  of  point  sources,  and  the  "mapping"  of  complex  events. 

It  is  the  purpose  of  this  contribution  to  give  results  for  the  excitation  of  seismic  waves  by  sources 
possessing  non-negligible  glut  moments  which,  it  is  hoped,  will  enable  certain  properties  of 
these  second-moment  tensors  to  be  estimated  from  the  data.  This  development  is  appealing, 
in  that  the  second-moment  tensors  contain,  in  a  unified  way,  information  on  the  spatial  extent, 
duration,  and  propagation  velocity  of  a  seismic  source.  Backus  has  considered  in  detail  the 
properties  of  these  moment  tensors  which  are  theoretically  determinable  under  a  number  of 
different  sets  of  assumptions.  Our  numerical  experiments,  which  will  not  be  discussed  fully 
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here,  suggest  that  in  practice  there  are  other  indeterminacies  which  remain  to  be  investigated. 
Nevertheless,  we  feel  that  to  pose  questions  about  the  spatial  and  temporal  extent  of  a  seismic 
source  in  terms  of  its  second-moment  tensors  is  a  useful  exercise,  in  that  if  the  inner-product 
matrix  and  corresponding  right-hand  vector  are  once  calculated  for  the  estimation  of  the  second- 
moment  tensor  components  from  a  given  set  of  seismograms,  they  represent  (in  a  highly  con¬ 
densed  form)  the  sum  of  information  in  the  data  about  the  size  and  duration  properties  of  the 
event.  It  is  then  possible  to  use  the  calculated  inner-product  matrix  to  test  various  models  of 
or  assumptions  about  the  earthquake,  or  to  invert  the  data  under  different  assumptions. 

The  seismic  displacement  of  a  spherically  asymmetric  Earth  model  due  to  any  indigenous 
source  may  be  written 


u(x,  t)  =  £  s.(x)  C  dt'  f  d3x'h.(t— t')  e^k'(x')  :  r(x',  1 


(III- 20 ) 


(for  notation,  see  Ref.  21).  Expanding  h,  (t  —  t')  e  (x1)  about  a  fiducial  point  x'  =  x  in  the  source 

K  -  -  .24  8 

region  and  a  fiducial  time  t'  =  t  close  to  the  origin  time,  we  obtain 

u(x.  t,  =  V  sk(x)  {hk(t  -  ts)  [M<°*°>e'k>  ♦  *  £  “fe® ’*£ U’ 


+  hk(t-ts)-  j-  M(0«2)5  J*} 


(111-21) 


where  the  modal  strains  e.(k)  and  strain  gradients  e/,kLe/k*  are  those  evaluated  at  the  fiducial 

point.  The  superscripts  on  the  various  glut  moment  tensors  indicate,  respectively,  the  spatial 
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and  temporal  degrees  of  the  moments,  ’  e.g.. 


=  f  dt  C  d3xf  ..(x,t) 


-  r  dt  r 

upq  j_„  Jv 


d  X(Xp  ~  Xsp)  (Xq  “  Xsq>  rij(^'t) 


C1’  =£ dtL  d3x(t-v (xP- v  w*t] 


(111-22) 


Z  5  21 

and  others  similarly.  Following  Gilbert  and  Dziewonski  and  Dziewonski  et  ah,  Eq. (111-21) 
leads  to  a  straightforward  evaluation  of  synthetic  seismograms  in  terms  of  the  spherical  com¬ 
ponents  of  the  moment  tensors,  when  the  spherical  components  of  the  strains  and  strain  gradi¬ 
ents  are  known.  The  calculation  is  most  economical  if  performed  in  a  coordinate  system  in 
which  the  source  is  at  the  pole  (6  =  0).  The  relevant  formulas  for  e.^.e^  have  been  given  in 
Refs.  25  and  21.  It  remains  to  list  the  corresponding  formulas  for  the  second  derivatives  e{.j>pq, 
and  the  spherical  components  of  this  fourth-rank  tensor  are  listed  in  Table  III-l.  They  have 
been  obtained  using  the  formalism  of  Phinney  and  Burridge. 


1 


The  constants  k  are 
n 
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k  -  1  (iL±l  .  (i  +n);il/2 
Kn  2n  *•  4x  (i  —  n)l y 

The  quantities  s^  for  spheroidal  modes  and  t.  for  toroidal  modes  are  of  the  form 


(III-23) 


s.  =  a1U  +  a2r"1U  +  a3r'2U  +  a4r"3U 


+  a5V  +  a^V  +  a2r'2V  +  agr-3V 


ti  =  a3W  +  a2r_1W  +  a3r'2W  +  a4r~3W 


(III- 24) 


etc.,  where  U,  V,  and  W  are  scalar  eigenfunctions  evaluated  at  the  source  depth.  The  coeffi¬ 
cients  a^  are  listed  for  spheroidal  modes  in  Table  III-2  and  for  toroidal  modes  in  Table  III- 3. 
It  may  be  noted  that,  as  a  result  of  the  compatibility  relations. 


We  have 


e. . 


ij.ki  +  eki ,ij  '  6 ik.il  T  ejl,ik 


+  e. 


(Ill—  25) 


s17  +  S4  "  2s13 


t3  +  t5  -  2t12 


s3  +s5  =  2s12 


tl6  +  t?  =  ttg  +  tJ5 


s10  ‘  s20 


S16  +  S7  =  S18  +  S15 


In  Tables  III-2  and  III- 3 

L  =  1(1+  1) 

where  I  is  angular  order. 


J.  H.  Woodhouse 


D.  INTERMEDIATE-  AND  LONG- PERIOD  SOURCE  MECHANISMS 
OF  SEVERAL  RECENT  EARTHQUAKES 

22  21 
In  the  last  SATS,  we  described  how  the  approach  of  Dziewonski,  Chou,  and  Woodhouse 

can  be  extended  to  incorporate,  in  the  search  for  the  earthquake  source  mechanism,  informa- 

21 

tion  on  excitation  at  very  long  periods.  In  the  original  application,  only  the  intermediate- 
period  data  (~60  s)  on  excitation  of  the  minor  arc  body  waves  were  used.  Here,  we  wish  to  show 
solutions  for  two  major  earthquakes  (Mg  >  7)  and  also  for  a  cluster  of  events  in  southern  Greece, 
which  gives  an  opportunity  to  examine  our  location  procedure. 

One  of  the  important  features  of  the  centroid  method  is  that,  in  addition  to  yielding  a  moment- 
tensor  solution,  it  also  provides  information  on  the  epicentral  location,  depth,  and  origin  time 
of  the  best  point  source.  Because  we  view  the  earthquakes  through,  as  if,  a  somewhat  distorting 
lens  —  lateral  heterogeneities  —  the  apparent  position  of  the  source  may  not  coincide  with  its  true 
location.  Functionals  of  the  Earth's  structure  employed  to  locate  earthquakes  using  arrival 
times  of  P-waves  and  using  the  centroid  method  are  quite  different.  Therefore,  the  answers 
may  be  biased  in  a  different  way.  Because  the  centroid  method  is  primarily  used  to  obtain 
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TABLE  1 1 1-2 


SPHEROIDAL-MODE  PARAMETERS  USED  IN  TABLE  lll-l 


Mode 

U 

f'0 

r-2U 

r3U 

D 

— 

r'2V 

- 

r’3V 

n 

1 

Bfl 

1 

■ 

-2 

2 

-) 

2 

-2 

H 

■ 

m 

1 

-2 

-2 

4 

II 

1 

B 

IL  +  2 

— (L  +  2) 

-L 

2L 

S5 

1 

-2 

2 

S6 

1 

-2 

V 

-1 

2 

-1 

tl+> 

-L 

s8 

i 

s9 

-1 

-2 

iL  +  ) 

s10 

2 

— (L  +2) 

-L 

jL(L  +  6) 

sll 

i 

7 

-1 

1 

1 

2 

1 

2 

1 

-1 

SI2 

1 

2 

-1 

1 

2 

-2 

3 

S13 

-l 

jL  +2 

4 

— — L  -  2 

2  L 

* 

-L 

1 

1 

5 

SI4 

2 

2 

~2 

3 

iL+i 

4L  2 

1 

iL+2 

4  2 

5.  1 

S)5 

2 

2 

~  4  L  ~  2 

SJ6 

-2 

IL+Z 

4  2 

5.  3 

~4L  ”  2 

s17 

-l 

2 

-2 

L 

S1 8 

3 

~2 

3 

1 

“2 

-(L+D 

s19 

-2 

+I-L  +  2 

J20 

2 

-L-2 

jL(L  +  6) 
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TABLE  111-3 

TOROIDAL-MODE  PARAMETERS  USED  IN  TABLE  111-1 


Mode 

W 

r’w 

r*2W 

r_3W 

*2 

-i 

2 

-2 

*3 

-2 

4 

*5 

i 

-2 

2 

*6 

1 

-2 

*7 

-i 

IL+1 

-L 

*8 

1 

*9 

-2 

IL+1 

*10 

-L 

jL(L  +  2) 

*n 

1 

i 

1 

2 

2 

*12 

1 

2 

-2 

3 

*13 

*>• 

* 

i' 

*14 

1 

1 

5 

2 

*15 

1 

2 

J(L  +  6) 

_2  L  -  — 

4  2 

*16 

_iL  +  l 

4  2 

+iL-  — 

4l  2 

*18 

1 

2 

1 

-1 

*19 

-1 

2 

-1 


2 


Men 


i  ■ 

j  information  on  the  source  mechanism  and  the  answers  depend  on  the  assumed  position  of  the 

j  source  (see  synthetic  example  in  Ref.  21),  it  is  possible  to  think  of  the  role  of  our  relocation 

procedure  as  the  means  to  minimize  in  the  least- square  sense  the  first  moment  of  the  stress 
23  24 

glut.  *  This,  regardless  of  the  physical  meaning  of  the  spatio-temporal  shift  in  the  apparent 
position  of  the  source,  leads  to  better  estimates  of  the  zeroth  moment,  commonly  referred  to 
as  the  seismic  moment  tensor.  Without  providing  a  rigorous  proof  of  this  statement,  we  shall 
show  an  example  that  illustrates  its  apparent  validity. 

1.  The  El  Asnam,  Algeria  Earthquake  of  10  October  1980 

This  event,  with  a  surface- wave  magnitude  of  7.3,  caused  the  death  of  at  least  5000  people. 
The  reported  length  of  the  surface  rupture  was  42  km.  Our  solution,  obtained  using  both  GDSN 
(Global  Digital  Seismograph  Network)  and  IDA  (International  Deployment  of  Accelerometers) 
data,  is  summarized  in  Table  III-4  and  the  nodal  plane  solution  corresponding  to  the  major 
double-couple  is  shown  in  Fig.  III-5.  The  seismic  moment,  for  the  focal  depth  of  10  km,  is 

about  4.5  x  lo2^  dyn-cm,  substantially  less  than  the  1.4  x  1027  reported  by  Baranowski  et  al.27 

28 

Although  our  moment  could  be  doubled,  if  a  sub-Moho  depth  were  to  be  assumed,  there  is  no 
reason  to  expect  that  much  of  the  radiation  was  due  to  the  fracture  at  sub-crustal  depths. 

Table  III-4  shows  the  moment-tensor  solution  obtained  using  the  starting  coordinates  (NEIS, 
10  s  added  to  allow  for  the  finite  duration)  and  the  final  centroid  location.  There  is  not  much 
difference  between  the  two,  because  the  position  of  the  epicenter  changed  very  little  (depth  was 
held  at  10  km)  —  less  than  20  km.  The  similarity  of  the  two  solutions  can  also  be  appreciated 
in  Fig.  Ill- 5. 

Figures  III-6(a)  through  (c)  and  III-7  show  comparison  of  the  observed  and  synthetic  seis¬ 
mograms  for  the  GDSN  and  IDA  data,  respectively.  The  match  for  the  GDSN  data  is  highly 
satisfactory,  and  it  is  difficult  to  imagine  that  our  estimate  of  the  moment  could  be  grossly 
wrong.  In  Fig.  III-7,  the  match  for  station  ERM  (Japan)  is  poor,  but  this  station  is  close  to 
the  nodal  plane  for  Rayleigh  waves.  When  the  observed  amplitudes  are  large  (TWO,  RAR, 

CMO),  the  agreement  between  the  observed  and  computed  seismograms  is  highly  satisfactory. 

2.  The  Eureka,  California  Earthquake  of  8  November  1980 

This  event,  given  an  Mg  of  7.2  by  the  NEIS,  occurred  offshore  and  caused  very  little  dam¬ 
age,  considering  its  magnitude.  It  is  probably  the  largest  earthquake  in  California  since  the 
Kern  County  series  of  19  52. 

Figure  III-8  shows  the  distribution  of  the  stations  used  to  study  this  event.  Table  III- 5 

gives  the  numerical  details  of  the  solution,  and  Fig.  III-9  shows  the  nodal  plane  solutions  for 

the  major  double-couple  for  the  starting  (NEIS)  and  final  (centroid)  locations  of  the  source. 

Unlike  for  the  El  Asnam  earthquake,  the  difference  between  the  two  solutions  is  significant. 

The  position  of  the  centroid  is  shifted  nearly  80  km  due  West,  and  this  change  in  location  has 

significant  impact  on  the  source  geometry  derived  in  inversion.  The  final  solution  is  a  nearly 

perfect  strike-slip,  in  accordance  with  the  tectonic  evidence.  The  starting  solution,  on  the 

other  hand,  has  a  very  large  dip-slip  component.  Our  final  solution  agrees  very  closely  with 

29 

that  obtained  for  this  earthquake  by  Lay  et  al. 

It  is  quite  clear  that  the  NEIS  location  was  not  wrong  by  80  km.  Also,  there  is  an  indication 
from  the  distribution  of  aftershocks  that  the  fault  was  bilateral.  Thus,  the  apparent  change  in 
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Principal  Axes 

P-Axis  N-Axis  T-Axis 


TABLE  III— 4 

ALGERIAN  EARTHQUAKE  OF  10  OCTOBER  1980t 


Starting 

Final 

Origin  Time 

12:25:23.5+10 

12:25:33 .5 

Latitude  (#N) 

36.20 

36.20 

Longitude  (°E) 

1.35 

1.58 

Depth  (km) 

10 

00) 

M 

3.74  ±  0. 08 

3.41  ±  0.05 

rr 

Mee 

-2.65  ±0.13 

-2.48  ±  0.09 

M 

-1.09±  0.21 

-0. 8 5  ±0.07 

Mre 

-1 . 10  ±  0.21 

-0.85  ±0.17 

M 

4 

-1.16  ±0.20 

—  1 . 58  ±0,17 

M9* 

-2.53  ±0.08 

-2.51  ±0.06 

Moment 

4.24 

3.92 

Plunge/AZM 

74°/102° 

720/9,0 

Moment 

0.68 

0.74 

Plunge/AZM 

10®/230° 

1 30/2290 

Moment 

-4.92 

-4.67 

Plunge/AZM 

120/322® 

ll®/322° 

26 

f  Scale  factor  is  10  dyn-cm. 


TABLE  111-5 

EUREKA  EARTHQUAKE  OF  8  NOVEMBER  1 980+ 


Starting 

Final 

Origin  Time 

1 0:27:33. 4H0 

10:27:45.6 

Latitude  (°N) 

41.16 

41.21 

Longitude  (°W) 

124.34 

125.40 

Depth  (km) 

14 

(14) 

M 

rr 

2.78  ±  0. 14 

-1.27  ±  0. 14 

Mee 

-5.90  ±  0.21 

-8.61  ±  0.20 

M 

99 

5.62  ±  0.19 

9.87  ±  0. 18 

Mre 

5.46  ±  0.28 

0.22  ±  0.38 

M 

r* 

6.97  ±0.28 

0.19  ±  0.46 

% 

2.58  ±0.16 

4.30  ±  0. 17 

T  -Ax  is 

Moment 

Plunge/AZM 

11.99 

36°/290° 

10.84 

,o/282= 

5 

< 

Moment 

-2.76 

-1.27 

i 

u 

c 

Z 

Plunge/AZM 

36°/52° 

88°/55° 

£ 

«/> 

X 

Moment 

-9.23 

-9.57 

r 

Q- 

Plunge/AZM 

34°/172° 

1  •/!  92° 

_ 

26 

t  Scale  factor  is  10  dyn-cm. 
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position  is  an  artifact  of  lateral  heterogeneity.  Yet,  shifting  the  epicenter  to  this  location  helps 

derive  the  correct  moment  tensor.  In  particular,  the  anomalously  large  elements  M  .  and 

r0  ru 

become  zero,  for  all  practical  purposes.  This  is  important,  because  these  two  elements  yield 
notoriously  unstable  values  for  shallow  earthquakes.3^  It  has  been  thought  that  the  small  values 
of  the  tangential  strain  near  the  surface  are  the  cause  of  this  instability.  In  a  way  this  is  true, 
because  the  tangential  stress  must  vanish  at  the  free  surface.  But  our  result  indicates  that 
lateral  heterogeneity  seriously  compounds  the  problem.  The  explanation  of  this  effect  follows. 

In  the  epicentral  coordinate  system,  excitation  of  seismic  waves  due  to  the  elements  M  „ 

1  4 

and  (dip-slip  components)  is  governed  by  the  spherical  harmonic  or  (9/90)  Y^  .  In  the 

asymptotic  limit,  these  are  shifted  in  phase  by  90°  with  respect  to  Y®,  Y^  or  their  derivatives, 
which  enter  into  the  expressions  for  excitation  due  to  the  remaining  four  independent  components 
of  the  moment  tensor.  Thus,  introduction  of  a  finite  (dip-slip)  fault  motion  may  mimic  a  phase 
shift.  Because  the  amplitude  due  to  these  terms  varies  as  cos  <p  or  sin  <p,  the  amount  of  the 
phase  shift  will  vary  like  cos  (0  +  e).  This  is  similar  to  the  phase  shift  introduced  by  a  change 
in  the  epicentral  coordinates  of  the  point  source,  or  the  equivalent  velocity  anomaly.  Yet, 
the  latter  effect  does  not  produce  the  amplitude  effects  that  introduction  of  and  Mr^  does. 
Therefore,  in  inversion  which  allows  for  shift  in  locations  of  the  source,  the  effect  of  a  true 
or  apparent  shift  of  the  epicenter  can  be  separated  from  the  excitation  pattern  due  to  the  dip- 
slip  elements  of  the  moment  tensor.  On  the  other  hand,  in  inversion  in  which  epicentral  co¬ 
ordinates  are  fixed,  as  they  are  in  all  other  inversion  techniques  and  in  the  starting  iteration 
of  our  method,  the  effect  of  lateral  heterogeneity  may  be  readily  absorbed  by  an  artificial  con¬ 
tribution  to  Mr0  and  M^. 

The  large  apparent  westward  shift  of  the  epicenter  of  the  Eureka  earthquake  is  not  an  ex¬ 
ception;  a  nearly  identical  shift  is  obtained  for  the  Mexicale  earthquake  of  1 5  October  1979. 

It  may  be  possible,  therefore,  to  anticipate  this  correction  for  the  strike- slip  earthquakes  in 
this  region  of  the  west  coast  of  North  America. 

Figures  IH-lO(a)  through  (c)  and  III  — 4 1  show  comparison  between  the  observed  and  synthetic 
seismograms  for  the  Eureka  earthquake. 

3.  Three  Earthquakes  in  Southern  Greece 

An  earthquake  of  magnitude  Mg  =  6.7  on  24  February  1981  was  followed  by  two  other  large 
shocks:  Mg  =  6.4  on  25  February,  and  6.5  on  4  March.  The  PDE  locations  indicate  that  all 
these  events  are  within  20  to  30  km  of  each  other,  and  it  occurred  to  us  that  this  may  represent 
an  opportunity  to  test  the  relative  accuracy  of  our  relocation  procedure.  Because  these  events 
were  so  closely  spaced  in  time,  the  network  coverage  was  essentially  constant  for  all  three 
earthquakes.  Only  the  GDSN  network  was  used;  IDA  data  are  not  yet  available.  Table  II 1  - 6 
lists  the  results  of  inversion  and  Fig.  Ill- 1 2  shows  the  fault-plane  solution,  indicating  that  all 
three  events  were  associated  with  normal  faulting. 

Figure  111-13  is  a  comparison  of  PDE  locations  (crosses)  with  the  relocation  result  (solid 
circles).  There  is  an  overall  shift  of  locations  toward  the  southwest;  this  may  well  be  a  large- 

scale  effect  of  the  lateral  heterogeneity:  the  Yugoslavian  earthquake  of  15  April  1979  and  the 

22 

Italian  event  of  23  November  1980  showed  westward  and  southwestward  shifts,  respectively. 

In  addition,  however,  the  relative  positions  of  individual  events  have  changed  after  inversion. 
While  the  PDE  locations  do  not  show  any  consistent  pattern  with  respect  to  the  strike  of  the 
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TABLE 

111-6 

RELOCATION  PARAMETERS  AND  THE  ELEMENTS 

OF  THE  MOMENT  TENSOR 

FOR 

THREE  RECENT  EARTHQUAKES  IN  SOUTHERN  GREECE' 

24  February  1981 

25  February  1  981 

4  March  1 981 

Origin  Time 

20:53:38.7 

2:35:53.3 

21:58:07.5 

LU  , 

7 

Latitude  (°N) 

38.24 

38.11 

38.31 

Longitude  (°E) 

23.02 

23.11 

23.18 

c 

.0 

% 

7.9 

5.4 

6.3 

8 

Latitude  (°N) 

37.91 

37.93 

38.01 

4) 

Longitude  (°E) 

22.60 

23.00 

22.98 

M 

rr 

-1.21  ±0.03 

-0.39  ±0.01 

-0.33  ±0.01 

M90 

1.40  ±  0.04 

0.40  ±0.01 

0.32  ±0.02 

M 

-0.20  ±0.03 

-0.01  ±  0.01 

0.01  ±  0.01 

99 

Mr0 

0.00  ±  0.06 

-0.06  ±0.02 

0.00  ±0.03 

M 

-0.45  ±0.04 

-0.05  ±0.02 

-0.11  ±0.02 

r9 

-0.06  ±0.03 

0.16  ±0.01 

0.12  ±0.01 

2 

Moment 

1.41 

0.46 

0.36 

< 

t— 

Plunge/AZM 

O 

* 

O 

5°/160o 

3°/160° 

Q> 

4! 

v> 

Moment 

-0.03 

-0.07 

0.01 

8. 

c 

z 

Plunge/AZM 

21  °/93° 

4°/70° 

17°/69° 

£ 

V) 

Moment 

-1.38 

-0.39 

-0.37 

< 

1 

CL. 

Plunge/AZM 

69°/271 0 

84°/301° 

73°/261° 

tThe  hypocentral  depth  has  been  held  fixed  at  15  km. 
is  common  for  all  moment  values. 

The  sea 

26 

e  factor  of  10  dyn-cm 
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fault,  the  events  of  25  February  and  4  March  are  located  along  the  strike  axis  of  the  initial  earth¬ 
quake.  Unless  this  is  a  coincidence,  the  conclusion  ^  aould  be  that  one  may  obtain  quite  precise 
relative  locations  by  matching  the  waveforms  of  GDSN  recordings;  data  from  fourteen  SRO, 

ASRO,  and  DWWSSN  stations  were  used  in  inversion. 

A.  M.  Dziewonski 


E.  PARAMETER  SPACE  SEARCH  FOR  THE  LOCATION 

OF  THE  BEST  POINT  SOURCE 

21 

Dziewonski,  Chou,  and  Woodhouse  have  developed  a  procedure  that  allows  simultaneous 
determination  of  the  moment  tensor  and  estimation  of  the  location  of  the  best  point  source  using 
long-period  (~ 60  s)  waveform  data  from  the  GDSN.  Because  their  method,  like  all  other  loca¬ 
tion  algorithms,  involves  linearization  of  the  problem,  the  convergence  of  the  solution  depends 
on  the  condition  that  the  starting  hypocentral  coordinates  be  linearly  close  to  the  true  position 
of  the  source  in  space  and  time.  After  having  obtained  satisfactory  results  for  numerous  sources 
(see  Ref.  21),  we  have  encountered  an  earthquake  for  which  it  seemed  impossible,  starting  with 
the  NEIS  hypocentral  coordinates,  to  obtain  a  solution  that  would  give  a  good  fit  to  the  data  and 
yield  a  moment  tensor  consistent  with  the  geological  and  geophysical  evidence.  This  was  the 
St.  Elias  earthquake  of  28  February  1979. 

The  inference  was  that  this  must  have  been  a  complex  event  of  considerable  spatio-temporal 
dimensions.  As  the  hypocentral  coordinates  determined  from  the  times  of  the  P-wave  arrivals 
led  to  unrealistic  answers,  the  alternative  was  to  perform  a  parameter  space  search.  Such  a 
search,  if  the  procedure  of  Ref.  21  were  to  be  applied  to  each  starting  solution,  could  be  ex¬ 
tremely  time  consuming.  The  subject  here  is  an  algorithm  that  allows  us  to  speed  up  this  pro¬ 
cess  by  a  factor  of  100  or  more. 

Given  the  origin  time,  epicentral  coordinates,  and  depth,  the  inverse  problem  for  the  mo¬ 
ment  tensor  of  a  point  source  is  derived  from  the  following  equations  of  condition; 


6 

uk(r,  t)  =  £  ^ki(I’rs*t)  '  fi  (III-26) 

i=l 

where  u^  is  the  record  in  the  set  of  seismograms  used  in  inversion,  are  the  excitation 
kernels,  and  f.  are  the  elements  of  the  moment  tensor. 

The  least-squares  condition  leads  to  normal  equations  of  the  following  form: 

A  ■  f  =  b  (lit—  27 ) 

where  an  element  of  the  inner  product  matrix  A  is: 


Aii 


*ki(t)  *kj(t)  dt 


(III-28) 


and  an  element  of  the  b  vector  is: 

U 


bj  =  2  [  k  «k(t) '  Vt}  dt 

i.  * i 


k  ll 


(IIT-29) 
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There  is  an  important  difference  between  the  changes  of  the  matrix  A  and  vector  b  with 
respect  to  perturbations  in  the  location  and  the  origin  time:  the  matrix  A  change  owly,  or 
is  relatively  insensitive  to  moderate  changes  in  hypocentral  parameters;  the  vector  b,  on  the 
other  hand,  changes  rapidly. 

The  general  proof  of  the  statement  above  may  be  difficult,  but  the  following  example  seems 
convincing.  Assume  that  we  consider  an  isolated  P-wave  with  the  accompanying  pP  and  sP 
pulses  with  a  dominant  period  of  1  s.  Let  us  say  that  our  station  is  at  a  distance  80'  from  the 
epicenter.  If  we  now  change  the  distance  to  81°,  the  integrals  (III- 2 8 )  will  change  very  little, 
because  the  synthetic  pulses  will  be  shifted  by  the  same  amount  of  time,  and  amplitudes  of 
P-waves  at  80'  change  very  slowly  with  distance.  At  the  same  time,  there  will  be  a  very  sig¬ 
nificant  change  in  the  result  of  calculation  of  the  elements  of  the  b  vector,  as  the  observed  and 
synthetic  pulses  will  be  shifted  with  respect  to  each  other  by  more  than  5  s.  It  is  possible, 
therefore,  to  think  of  a  scheme  in  which  we  could  calculate  the  matrix  A  only  once  and  assume 
that  its  value  remains  constant  with  respect  to  moderate  changes  in  location  of  the  source,  but 
b  would  have  to  be  calculated  exactly. 

Let  us  consider  what  enters  into  the  calculation  of  b.  The  excitation  functions  if  can  be 
written  [see  Eqs.  (A15)  through  ( A 1 7 )  in  Ref.  21]: 


‘ik  =  E  E  <Tj(,)<t.rs>  ^‘(0,0.0  (III- 30) 

•  3 

where  r  .  is  the  source  radius;  the  quantities  a  represent  contributions  of  the  appropriate  com¬ 
binations  of  eigenfunctions  and  their  derivatives  for  all  overtones  of  a  particular  angular  order 
number  l . 

The  vector  b  is  therefore: 


Vi  ES 

t  k  J 


where 


rir-i,  k  V '/'W-1  " 


(III-  31 ) 


(III— 32) 


The  important  property  of  expression  (II I  -  32 )  is  that  it  does  not  depend  on  t*  e  epicentral 
coordinates.  It  is  possible,  therefore,  to  calculate  y.(J*  only  once  for  a  giv  i  depth  and  origin 
time  and  then  evaluate  Eq.  (Ill- 31)  very  rapidly  for  anv  position  of  the  source,  as  only  a  sum  of 
scalar  values  is  involved,  while  in  evaluation  of  Eq .  (Ill-  12)  one  must  consider  100  to  200  sam¬ 
ples,  on  average.  Such  a  transformation  is  not  possible  in  the  case  of  matrix  A.  which  is  qua¬ 
dratic  in  functions  that  depend  on  the  position  of  the  source. 

The  procedure  is  as  follows.  We  establish  a  grid  system  in  the  geographical  coordinates, 
calculate  y  V  *  and  A  for  the  center  of  the  grid  (r  )  and  then,  very  rapidly,  find  b(r)  for  each 

JK  —  — C  ~  ~ 

grid  point.  The  solution  is  then  found: 


f(r)  b(r)  •  4_1(rc) 


(111  —  33) 
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The  goodness  of  fit  is  approximated  by  the  relative  variance: 


e 


where 


(III- 34) 


It  is  important  to  remember  that  the  vector  f(r)  represents  only  an  approximate  solution; 
thus,  it  is  possible  that  this  estimate  of  the  variance  may  become  negative. 

Figure  III—  1 4  shows  a  result  of  application  of  the  method  to  a  synthetic  case,  in  which  a 

point  source  (normal  fault)  has  been  placed  at  the  center  of  the  grid;  the  grid  spacing  is  0.2° 

of  arc.  The  network  coverage  is  realistic:  the  same  as  that  used  in  the  investigation  of  the 
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trpinia  earthquake  of  23  November  1980.  *  The  value  of  e  (the  figures  indicate  100  x  e)  in¬ 

deed  reaches  a  minimum  at  the  location  of  the  point  source  (cross).  The  noticeable  distortion 
of  contours  is  the  result  of  uneven  network  coverage. 

In  the  second  synthetic  example,  we  move  the  source  by  1  “of  longitude  and  i'  of  latitude 
from  the  center  of  the  grid.  Figure  III— 1 5  shows  the  result.  The  center  of  the  minimum  value 
of  e  is  only  a  fraction  of  the  grid  spacing  away  from  the  "true"  position  (cross),  even  though 
this  point  is  some  150  km  away  from  the  grid  center  (plus  sign).  However,  e  assumes  negative 
values;  this  is  the  result  of  the  approximation  implied  by  Eq.  (Ill- 33).  This  means  that  the  ele¬ 
ments  of  the  moment  tensor  are  somewhat  incorrect.  But  the  important  result  is  that  the  loca¬ 
tion  of  the  source  has  been  found  correctly.  Having  this  location,  we  can  initiate  the  procedure 
of  Ref.  21,  and  the  exact  result  will  be  obtained  in  one  or  two  iterations. 

Shifting  the  source  in  time  may  also  produce  rather  dramatic  effects.  Because  we  deal  with 
signals  in  which  energy  is  concentrated  within  the  range  of  periods  from  50  to  60  s,  shifting  the 
origin  time  by  25  to  30  s  may  lead  to  an  apparent  convergence  of  the  solution  but  with  a  wrong 
sign  of  the  moment  tensor;  these  secondary  minima  are  quite  well  defined. 

Thus,  for  complex  earthquakes  of  significant  duration  we  produce  a  series  of  contour  maps, 
each  for  different  offset  in  the  origin  time;  a  step  between  5  and  10  s  is  reasonable  for  the  fre¬ 
quency  band  of  our  analysis.  Application  of  this  procedure  to  the  recordings  of  the  St.  Elias 
earthquake  by  the  SRO/ASRO  network  allows  us  to  locate  the  position  of  the  best  point  source 
some  100  km  east  of  the  epicenter  and  37.5  s  after  the  origin  time  determined  from  the  first 
breaks  of  the  P-waves.  This  is  shown  in  Fig.III-16. 

Having  determined  the  position  and  time  for  the  best  point  source,  we  can  now  apply  the 
22 

procedure  of  Dziewonski  for  the  analysis  of  multiple  sources.  Since  the  St.  Elias  earthquake 
generated  P-waves  of  significant  amplitude,  it  is  reasonable  to  assume  for  the  first  source 
the  origin  time  and  location  of  die  NEIS  epicenter,  and  for  the  second  source  the  result  shown 
in  Fig.III-16.  These  times  and  locations  change  surprisingly  little  during  the  inversion. 
Figure  III-17  is  a  map  taken  from  Stephens  et  al.32  The  fault-plane  solutions  corresponding  to 
the  "major  double- couple"  derived  from  the  analysis  of  the  moment  tensor  are  shown  at  the 
bottom  of  the  figure. 
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The  first  source  is  a  very  shallow  thrust,  in  agreement  with  the  nodal-plane  solution  of 
33 

Perez  and  Jacob  based  on  the  signs  of  the  P-wave  arrivals.  The  steeply  dipping  nodal  plane 
of  these  authors  coincides  with  that  of  Fig.  111—17  within  a  few  degrees;  the  dip  direction  of  the 
shallow  dipping  plane  is  usually  poorly  determined.  But  the  second  source  with  the  centroid 
time  of  37  s  after  the  origin  has  a  substantially  different  mechanism.  It  is  a  reverse  fault, 
with  the  fault  planes  dipping  at  roughly  45°  and  a  substantial  strike- slip  component.  The  sense 
of  slip  is  consistent  with  the  north- south  compression.  Thus,  our  result  implies  that  the  source 
mechanism  of  the  St.  Elias  earthquake  underwent  a  major  change  during  the  faulting  process. 

Using  our  procedure,  we  have  also  analyzed  the  Tabas-e-Golshan  (Iran)  earthquake  of 
16  September  1978.  The  best  point  source  was  located  northeast  of  the  epicenter  and  some 
20  s  later.  Unfortunately,  the  network  coverage  for  this  event  is  very  uneven,  and  there  are 
major  trade-offs  between  the  geographical  coordinates  and  the  origin  time.  We  have  fixed  the 
locations  of  times  of  three  sources,  as  shown  in  Fig.  III-18,  and  simply  solved  for  their  source 
mechanism.  Again,  for  the  first  source  we  obtain  an  essentially  dip-slip  solution,  while  the 
second  and  third  sources  are  reverse  faults.  As  in  the  case  of  the  St.  Elias  event,  a  significant 

strike- slip  component  develops  near  the  end  of  the  faulting  process.  The  initial  solution  is 
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essentially  consistent  with  that  obtained  by  Berberian  et  al.  from  the  P-wave  polarities. 

The  techniques  described  in  this  report  provide  a  tool  for  resolving  very  rapidly  the  es¬ 
sential  complexities  of  large  earthquakes.  If  the  data  were  transmitted  in  real  time,  the  results 
for  the  February  1979  or  September  1978  events  could  be  obtained  within  3  h  after  the  earth¬ 
quakes  occurred.  This  could  be  important  for  rapid  assessment  of  the  damage,  as  well  as  for 
determining  the  strategy  of  post-seismic  geophysical  and  geodetic  measurements. 

A.  M.  Dziewonski 
J.  H.  Woodhouse 


F.  PURE- PATH  MODELS  FROM  RAYLEIGH- WAVE  DISPERSION 

IN  A  FREQUENCY  RANGE  FROM  1.5  TO  6  mHz 

35  3  £> 

Dziewonski  and  Dziewonski  and  Steim  have  described  a  new  technique  of  analysis  of 
mantle  waves  which  allows  one  to  extract  directly  information  on  the  Earth's  structure  from 
waveform  data.  An  iterative  procedure  was  designed  to  perturb  the  structural  elements  (both 
velocity  and  Q)  to  improve  the  transfer  function  for  Rayleigh  waves  traveling  around  the  world. 
We  have  presented  examples  of  application  of  this  technique  in  the  analysis  of  synthetic  as  well 
as  real  data.  We  have  named  this  method  the  "waveform  inversion  technique  (WIT)." 

This  report  describes  application  of  the  method  to  a  "pure-path"  decomposition  of  the  phase- 
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velocity  and  Q  data.  Numerous  analyses  of  this  type  have  been  performed  in  the  past. 

The  results  obtained  in  these  studies  have  been  used  to  infer  the  depth  to  which  lateral  varia¬ 
tions  in  the  elastic  and  anelastic  parameters  are  correlated  with  the  surficial  tectonic  nature 
of  a  given  region.  Implications  of  these  inferences  can  be  of  fundamental  importance  to  our 
understanding  of  the  dynamics  of  the  Earth.  The  results  presented  in  this  section  might  be 
understood  as  a  test  of  WIT  in  application  to  this  type  of  analysis.  Also,  there  is  an  indica¬ 
tion  that  extension  of  measurements  to  frequencies  as  low  as  1.5  mHz  provides  important 
structural  constraints. 
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From  a  large  number  of  visually  examined  recordings  of  the  IDA  network,  we  have  se¬ 
lected  37  recordings  that  showed  the  best  SNR.  Their  northern  hemisphere  poles  are  shown 
in  Fig.  III-19.  The  poles  have  a  fairly  uniform  distribution  in  longitude,  but  nearly  equatorial 
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latitudes  (corresponding  to  nearly  polar  paths)  are  poorly  represented.  In  the  future,  particular 

effort  should  be  made  to  assure  even  geographical  coverage. 

The  matching  of  transfer  functions  for  individual  paths  was  achieved  by  perturbation  of  three 

parameters  in  shear  velocity  (constant  perturbations  in  the  depth  ranges  80  to  220,  220  to  400, 

and  400  to  670  km)  and  two  in  Q  (perturbations  in  the  depth  ranges  from  80  to  400  and  400  to 
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670  km).  The  experiments,  such  as  those  discussed  in  our  previous  reports/  '  demonstrated 
that  this  parametrization  is  sufficient  to  obtain  good  transfer  functions.  Clearly,  other  parame- 
trization  schemes,  including  transverse  anisotropy are  feasible.  The  range  of  frequencies 
(angular  order  numbers)  used  in  inversion  was  determined  through  visual  examination  of  the 
spectra:  the  total  span  of  the  order  numbers  used  is  from  7  to  60  (from  800  to  150  s,  roughly), 
but  it  varied  appreciably  from  record  to  record.  For  this  reason,  one  should  not  anticipate 
full  correspondence  between  the  results  of  regionalization  using  either  the  model  or  data  space. 

The  regionalization  adopted  here  is  based  on  the  work  of  Mauk  who  divided  the  world  into 
20  types  of  regions,  using  a  5°  x  5“  grid,  and  considered  fractional  contributions  (with  a  pre¬ 
cision  of  one-tenth)  of  individual  regions  to  a  given  cell.  We  have  combined  Mauk's  regions 
into  four  types:  (1)  stable  continents,  including  continental  shelves,  for  which  tectonic  activity 
ceased  at  least  400  m.y.  ago;  (2)  tectonic,  including  island  arcs;  (3)  ocean  floors  younger  than 
38  m.y.;  and  (4)  ocean  floors  older  than  38  m.y.  Perhaps  a  more  detailed  regionalization  will 
be  warranted  in  the  future,  when  the  number  of  paths  increases  significantly.  Figure  III— 20  is 
a  map  of  the  world  schematically  showing  our  regionalization.  Fractional  contributions  for 
each  cell  were  considered  in  actual  calculation  of  path  compositions. 

After  the  data  obtained  through  application  of  WIT  to  individual  paths  were  corrected  for 
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ellipticity,  *  the  pure-path  decomposition  was  calculated  using  the  approach  of  Toksoz  and 
37 

Anderson.  The  results  for  periods  of  normal  modes,  group  velocities,  and  Q  are  presented 
in  Table  III-7  for  the  range  of  angular  order  numbers  from  9  to  55.  The  differences  of  pure- 
path  phase  and  group  velocities  from  those  computed  for  PREM1^  are  shown  in  Fig.  III-2f(a-b). 
Roughness  of  the  curves  at  frequencies  lower  than  3  mHz  and,  in  particular,  higher  than  5  rrTTz 
is  the  result  of  the  fact  that  data  were  not  available  for  all  paths  outside  the  3-  to  5-mHz  range. 
Except  for  roughness,  there  is  no  overall  change  in  the  trend  due  to  variations  in  path  coverage; 
this  indicates  the  stability  of  the  inverse  problem  and  internal  consistency  of  the  data  set. 

At  frequencies  higher  than  4  mHz,  the  major  differences  are  between  young  oceans  and  tec¬ 
tonic  regions  on  one  hand  and  stable  continents  and  old  oceans  on  the  other.  This  is  qualita¬ 
tively  consistent  with  the  results  of  previous  studies  of  this  kind.  Phase  velocities  become  very 
close  to  each  other  at  the  lower  frequency  range  of  analysis  (1.57  mHz),  but  then  begin  to  diverge 
in  a  rather  distinct  manner:  for  stable  continents  and  tectonic  regions  the  phase  velocities  follow 
the  global  average  and  are  very  close  to  each  other  up  to  a  3-mHz  frequency,  while  the  values 

for  young  oceans  begin  to  diverge  much  earlier.  A  similar  trend  can  be  inferred  from  the  re- 
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suits  of  Silver  and  Jordan  for  the  end  members  of  their  regionalization  of  continental  and 
oceanic  areas.  An  intuitive  interpretation  would  imply  deep-seated  differences  in  structure 
beneath  the  young  and  old  oceans. 

The  results  for  pure-path  Q  indicate  that  the  derived  differences  are  not  statistically  sig¬ 
nificant;  there  is  some  indication,  however,  that  Q  for  the  continental  regions  may  be  higher 
than  for  the  oceanic  ones.  We  test  this  hypothesis  by  decreasing  the  resolution  of  the  regionali¬ 
zation  and  distinguish  only  between  continents  and  oceans.  The  results  are  shown  in  Table  III-8 
which  also  contains  the  properly  weighted  global  average  and  comparison  with  the  appropriate 
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TABLE  1 11—7 

RESULTS  OF  PURE-PATH  ANALYSIS  FOR  PERIODS  OF  NORMAL  MODES 
GROUP  VELOCITIES,  AND  ATTENUATION 
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functionals  for  PREM.15  The  standard  errors  in  pure-path  Q  decrease  substantially  and  could 
be  considered  statistically  significant.  The  systematic  trend  of  differences  between  the  conti¬ 
nental  and  oceanic  Q  values  is  illustrated  in  Fig.  III-22.  At  a  period  of  250  s,  the  continental 
and  oceanic  Q's  are  196  ±  16  and  160  ±  7,  respectively  —  a  difference  of  about  20  percent. 

The  comparison  of  the  average  data  from  37  paths  studied  in  this  paper  with  the  functionals 
of  PREM,  a  model  derived  using  a  much  broader  data  set,  is  a  test  of  the  overall  performance 
of  WIT.  Even  small  systematic  differences  could  be  detected  in  this  way.  For  periods  of  nor¬ 
mal  modes,  the  maximum  difference  is  0.05  percent  with  an  average  rms  of  0.02  percent  —  an 
excellent  agreement  indeed!  The  differences  between  the  group  velocities  do  not  significantly 
exceed  0.01  k m/s,  a  remarkable  result  considering  the  fact  that  group  velocities  were  not  used 
in  derivation  of  PREM  and  that  there  are  substantial  variations  in  group  velocities  estimated 
for  individual  paths:  the  maximum  difference  in  pure-path  group  velocities  approaches  5  percent. 
The  average  Q  values  are  also  very  close  to  those  for  PREM.  At  all  frequencies  the  difference 
is  less  than  5  percent,  which  is  the  estimated  accuracy  of  the  average  Q  data  in  the  relevant 
range  of  frequencies  (see  Table  V  in  Ref.  16). 

Because  of  the  inherent  difficulty  in  the  direct  comparison  of  the  results  of  measurements 
of  Q  with  the  original  data  (seismograms),  we  show  in  Fig.III-23  the  average  Q  values  obtained 
in  this  study  together  with  the  results  of  Sailor,50  and  Geller  and  Stein51;  the  Q-curve  for  PREM 
is  included.  Sailor  obtained  his  results  by  using  the  traveling-wave  approach  and  a  refined 
spectral  ratio  method;  thus,  there  should  be  no  common  source  of  bias  between  his  result  and 
ours.  Sailor's  data  show  small  amplitude  oscillations  superimposed  on  a  smooth  curve  that  is 
essentially  parallel  to  ours  with,  perhaps,  a  few-percent  systematic  difference  at  the  short- 
period  end.  The  results  of  Geller  and  Stein,  obtained  by  measuring  the  time-domain  decay  of 
the  envelope  of  individual  modes,  have  very  large  scatter;  this  might  be  construed  as  an  indica¬ 
tion  of  the  advantage  of  the  traveling-wave  approach  for  periods  shorter  than  about  500  s. 

Finally,  we  consider  the  results  obtained  in  the  model  space.  For  each  path  we  have  a  set 
of  perturbations  in  velocities  and  Q  that  were  introduced  to  improve  the  transfer  function.  We 
can  thus  solve  the  pure-path  problem  directly  in  the  model  space,  without  the  need  to  invert  the 
data.  As  stated  above,  the  results  need  not  be  fully  compatible  with  those  obtained  for  the 
normal-mode  functionals,  because  the  frequency  range  of  inversion  differs  from  path  to  path, 
while  the  parametrization  remains  the  same.  Figure  III- 24  shows  perturbations  in  shear  ve- 
locity  for  the  individual  regions  with  respect  to  the  reference  model  (IREM  ).  The  most  un¬ 
usual  feature  of  this  plot  is  the  large  differences  between  the  velocities  for  young  and  old  oceans 
in  a  depth  range  from  400  to  670  km.  At  shallower  depths,  the  results  are  compatible  with 
those  from  other  studies.  The  structures  in  the  depth  range  from  80  to  220  km  reflect  the 
composite  properties  of  the  crust,  lithosphere,  and  the  low-velocity  zone:  fundamental  Rayleigh- 
mode  data  with  a  short-period  cutoff  at  165  s  have  no  resolving  power  to  determine  these  sepa¬ 
rate  elements  of  the  structure. 

The  total  spread  of  values  is  very  large  (0.34  km/s)  in  the  depth  range  from  80  to  220  km. 
The  relative  slowness  of  the  young  oceans  and  tectonic  regions  is  in  good  agreement  with  re¬ 
gional  surface-wave  studies.  The  total  spread  of  shear  velocities  is  relatively  small  in  the 
depth  range  from  220  to  400  km.  It  is  less  than  0.1  km/s,  and  the  error  bars  in  Fig.  111-24 
indicate  that  the  differences  are  not  statistically  significant.  In  a  depth  range  from  400  to 
670  km,  the  velocities  for  stable  continents  and  tectonic  regions  are  very  close  to  each  other 
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and,  also,  to  the  global  average.  The  difference  between  the  young  and  old  oceans  Is  large  — 
nearly  0.2  km/s  or  4  percent  —  and  appears  to  be  statistically  significant.  One  cannot  discount 
the  possibility  that  it  is  an  artifact  of  the  errors  introduced  by  the  geometrical-ray-theory  ap¬ 
proximation,  but  in  that  case  abnormal  properties  of  the  models  for  stable  continents  and  tec¬ 
tonic  regions  could  be  expected:  either  of  the  two  region  types  is  smaller  in  the  total  area  than 
the  young  oceans. 

The  matter  should  be  satisfactorily  resolved  by  application  of  the  theory  of  Woodhouse  and 
r,irnius,1’i  which  avoids  the  approximation  inherent  in  the  geometrical  ray  theory.  Also,  im¬ 
proved  path  coverage  should  lead  to  higher  resolution.  The  important  result  of  this  experiment 
is  •bat  the  extended  period  range  provides  sufficient  resolution  below  a  depth  of  400  km  to  infer 
structural  differences  in  the  region  of  the  transition  zone. 

A.  M.  Dziewonski 
J.  M.  Steim ' 


t  Department  of  Geological  Sciences,  Harvard  University. 
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DISPLACEMENT  SPECTRUM 
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Fig.  Ill— 1.  Spectrum  of  P-wave  recorded  by 
short-period  SRO  station  ANMO  for  a  deep 
focus  earthquake  in  Banda  Sea.  Dotted  line 
is  least-squares  line  fit  to  spectrum  for  fre¬ 
quencies  greater  than  0.7  Hz.  Solid  lines  are 
predicted  spectral  decay  for  t£  =  0.4,  0.2,  and 
0. 1  for  an  equivalent  event  at  surface.  Arrows 
show  "time  acceptance  window"  of  P-wave. 


DISPLACEMENT  SPECTRUM 


Fig.  Ill— 2.  Spectrum  of  P-wave  recorded  by 
short-period  SRO  station  NWAO  for  a  deep 
focus  earthquake  in  Peru.  Lines  and  sym¬ 
bols  defined  in  caption  to  Fig.  Ill- 1 . 


FREQUENCY  (Hi) 


1 


Fig.  III-3.  Synthetic  ScS  and  ScP  phases 
computed  for  varying  values  of  the  lower 
time  constant  rm  of  a  mantle-wide  ab¬ 
sorption  band.  For  frequencies  less  than 
1/2  irrm,  the  attenuation  model  converges 
to  that  of  PREM;  for  frequencies  greater 
than  1/2  firm  Hz,  Q a  and  Qp  decrease  as 
to-1.  Results  shown  are  for  a  triangle 
source  pulse  of  1-s  width  convolved  with  a 
short-period  WWSSN  instrument  response. 
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Fig.  XII-4.  Synthetic  spectral  ratio  of  ScS/ScP 
vs  frequency.  Solid  lines  are  for  a  mantle 
wide  absorption  band  that  converges  to  PREN 
for  frequencies  less  than  1/2  irrm  Hz  and  hat 
Qq.  and  Q0  decrease  as  w'l  for  frequencies 
greater  than  1/2  irrm.  Dotted  lines  include  a 
low  Q  zone  in  the  lower  150  km  of  mantle  with 
Qq  =  100  and  Qp  =  50.  Dashed  lines  are  for 
an  absorption  band  spanning  the  band  of  high- 
frequency  body  waves  (tjh  =  0.01)  in  the  upper 
400  and  lower  150  km  of  the  mantle,  but  shifted 
to  low  frequencies  in  mid-mantle  lrm  =  0.6. 
0.4).  Shaded  region  is  domain  of  frequencies 
and  ScS/ScP  observations  found  by  Burdick/1 


Fig.  IH-5.  Comparison  of  fault-plane  solutions 
(major  double-couple)  obtained  using  starting  and 
final  coordinates  of  El  Asnam,  Algeria  earth¬ 
quake  of  10  October  1980. 
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Fig.  IU-8.  Coverage  of  Eureka  earthquake  of  8  November  1980 
by  instruments  of  GDSN  (SRO  -  square;  ASRO  -  upward -pointing 
triangle)  and  IDA  (diamonds)  networks. 
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Fig.  Ill -9.  Comparison  of  fault-plane  solu¬ 
tions  (major  double -couple)  obtained  using 
starting  and  final  coordinates  of  Eureka 
earthquake  of  8  November  1980. 
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Fig.  m-lO(a-c).  Comparison  of  observed  (top)  and  synthetic  GDSN 
seismograms  for  Eureka  earthquake  of  8  November  1980. 
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Fig.  m-lO(a-c).  Continued. 
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Fig.  Hl-lO(a-c),  Continued. 
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Fig.  111-11.  Comparison  of  observed  (top)  and  synthetic  IDA  seismograms 
for  Eureka  earthquake  of  8  November  1980. 


Fig.  111-12.  Fault-plane  solutions  (major 
double- couple)  of  three  closely  spaced 
earthquakes  in  Southern  Greece. 
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Fig.  III-13.  Relocation  of  epicenters  of 
three  earthquakes  in  Southern  Greece. 
Crosses  are  NEIS  locations,  solid  cir¬ 
cles  are  centroid  locations.  Note  that, 
after  relocation,  earthquakes  lie  on  an 
axis  parallel  to  strike  direction  (see 
Fig.  111-12). 
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Fig.  III-14.  A  contour  map  of  e  X  100  [Eq.  (III-34)].  Minimum  value  of  e 
indicates  best  location  of  point  source  for  specific  amount  of  shift  inorigin 
time.  In  this  test  with  synthetic  data,  location  of  source  coincided  with 
center  of  grid  system.  Result  is  exact. 
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Fig.  Ill—  1 5,  Same  as  Fig.  III-14,  but  source  is  offset  with  respect  to  center 
of  grid  (plus  sign).  Parameter  search  allows  locating  of  source  satisfactorily 
even  though  offset  is  greater  than  150  km. 
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Fig.  III-16.  Application  of  parameter  search  procedure  to  analysis  of 
St.  Elias  earthquake  of  28  February  1979.  Crossed  circle  is  NEIS  lo¬ 
cation  of  this  event.  We  determine  that  best  location  of  point  source 
is  about  100  km  to  east  of  NEIS  location  and  37  s  after  origin  time. 


Fig.  UI-17.  Result  of  inversion  of  SRO/ASRO 
from  St.  Elias  earthquake  for  point  sources: 
first  one  coincides  with  NEIS  location  (star), 
second  source  (cross)  corresponds  to  lo¬ 
cation  determined  in  Fig.  III-16.  Mechanism 
of  source  i  closely  corresponds  to  that  obtained 
from  P-wave  polarities^;  source  2  is  mark¬ 
edly  different.  Tectonic  map  is  from  a  paper 
by  Stephens  et  al.32  [Reprinted  with  permis¬ 
sion  from  C.  D.  "Stephens  et  al. ,  Bull.  Seismol. 
Soc.  Am.  70,  1607-1634  (1980).] 
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Fig.  m-18.  Result  of  inversion  of  SRO/ASRO  data  from  Tabas-e-Golshan 
(Iran)  earthquake  for  three  point  sources  denoted  by  crossed  circles  (se¬ 
quence  is  from  lower  right  to  upper  left).  Faulting  progresses  from  dip- 
slip  (or  very  shallow  thrust)  to  reverse  faulting,  and  last  source  develops 
significant  strike-slip  component.  Note  similarity  between  evolution  of 
this  event  and  that  of  St.  Elias  earthquake.  Map  is  from  Berberian  et  al.34. 
solution  from  P-wave  polarities  is  superimposed  on  map.  [Reprinted  with 
permission  from  M.  Berberian  et  aU,  Bull.  Seismol.  Soc.  Am.  69,  1851- 
1860  (1979).]  — 
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Fig.  in-19.  Location  of  northern  hemisphere  poles  of  37  great-circle 
paths  analyzed  in  this  report. 
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+  :  TECTONICALLY  ACTIVE  I:  OCEAN  FLOOR  YOUNOER  THAN  38  m.y. 

—  ■  STABLE  CONTINENTS  AND  SHELVES  BLANK:  OCEAN  FLOOR  OLDER  THAN  38m.y. 


Fig.  IE-20.  Schematic  representation  of  regionalization  adopted  in  our 
pure-path  analysis.  This  regionalization  has  been  obtained  by  mbining 
several  regions  of  similar  tectonic  nature  as  defined  by  Mauk.1 
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Fig.  Ill— 21.  Pure-path  phase  velocities  (a)  and  group  velocities  (b)  obtained 
through  decomposition  of  results  of  measurements  for  37  great-circle  paths. 
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Fig.  IU-22.  Attenuation  (Q)  for  continental  and  oceanic  regions;  data  are  not 
sufficiently  numerous  and  accurate  to  allow  four-region  analysis  for  Q.  Dif¬ 
ferences  between  continental  and  oceanic  Q  values  are  statistically  significant; 
see  Table  IIX-8. 
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Fig.  Ill— 23.  Comparison  of  Q  data  of  Sailor50  (traveling-wave  approach),  and 
Geller  and  Stein§f  (standing-wave  approach)  with  those  obtained  in  this  study; 
values  computed  for  PREM*°  are  added  for  reference. 
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GLOSSARY 


AA 

Automatic  Association 

ABAR 

Arrivals  Based  Analyst  Review 

ASL 

Albuquerque  Seismological  Laboratory 

ASRO 

Abbreviated  Seismic  Research  Observatory 

ATN 

Augmented  Transition  Network 

DARPA 

Defense  Advanced  Research  Projects  Agency 

DBS 

Database  Subsystem 

DEC 

Digital  Equipment  Corporation 

DMA 

Direct  Memory  Access 

DTW 

Dynamic  Time  Warping 

DWWSSN 

Digital  WWSSN  Stations 

FEL 

Final  Event  List 

GDSN 

Global  Digital  Seismograph  Network 

IDA 

International  Deployment  of  Accelerometers 

IESD 

International  Exchange  of  Seismic  Data 

IREM 

Interim  Reference  Earth  Model 

LTA 

Long  Term  Average 

NEIS 

National  Earthquake  Information  Service 

PDE 

Preliminary  Determination  of  Epicenters 

PEL 

Preliminary  Event  List 

PREM 

Preliminary  Reference  Earth  Model 

RAG 

Related  Arrival  Group 

RSTN 

Regional  Seismic  Test  Network 

SAS 

Seismic  Analysis  Station 

SATS 

Semiannual  Technical  Summary 

SDAC 

Seismic  Data  Analysis  Center 

SDC 

Seismic  Data  Center 

SNR 

Signal-to-Noise  Ratio 

SRO 

Seismic  Research  Observatory 

STA 

Short  Term  Average 

USGS 

U.S.  Geological  Survey 

UUCP 

UN  DC -UN  IX  Copy 

WIT 

Waveform  Inversion  Technique 

WWSSN 

World-Wide  Standard  Station  Network 
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